Optimize FIR filter kernel in Matlab

Question

I am currently using the firls function in Matlab to create my FIR filter kernel. This works well but I typically have to play with the filter order to get the gain to 1 and the shape to be somewhat Gaussian in the frequency domain.

I would like to write an algorithm that finds the best filter kernel given a range of passband frequencies for the FIR filter. My initial attempt was to create an ideal filter kernel, calculate the integral and then step through a range of orders and settle on the order which had an integral closest to the ideal kernel. The issue with this approach was the chosen filter orders were too high and the shape of the kernel was not a smooth gaussian as shown on the left.

Left graph shows desired kernel, right is undesired higher order kernel. Red point plot is the ideal filter. Y-axis is frequency.

Here is the Matlab code I have been using to generate the filter and what I'm hoping to expand off of to create a gaussian-like, smooth FIR filter for any inputted frequency range.

clc,clear 
close all


srate = 250e3;
nyquist = srate/2;
npnts = srate*5;    %generate number of points for 5 seconds of sampling
time = (0:npnts-1)/srate; 

% simulate brownian noise
bnoise = 30*cumsum(randn(1,npnts)); 
% simulate white noise 
wnoise = 50*randn(1,npnts);

% vector of frequencies in kHz
hz = linspace(0,srate/2,floor(npnts/2)+1)/1e3;  %go from 0 to nyquist, frequency resolution is defined by last term

% signal settings 
freq1 = 26e3;
ampl = 1;
Y = ampl*sin(2*pi*freq1*time) + bnoise + wnoise;



% amplitude spectrum via Fourier transform
signalX = fft(Y);
signalAmp = 2*abs(signalX)/npnts;   %need to multiply by 2 to recover amplitude from negative freqs
                                    % divide by npnts to normalize fourier coefficients 
                                   
% FIR filter specs
passband = [freq1-2e3 freq1+2e3];
transw = 0.01;   %how the filter's edges taper
test_order = 7;

order = round(test_order*srate/passband(1)); % define the number of time points for the filter kernel

shape = [ 0 0 1 1 0 0 ];    %FIR shape

%define frequency shape of the FIR shape, firls function requires
%frequencies to go from 0 to 1
frex = [0 passband(1)-passband(1)*transw passband passband(2)+passband(2)*transw nyquist]/nyquist;

%define kernel
filter_kernel = firls(order, frex, shape);

% power spectrum of filter kernel
filtkernX = abs(fft(filter_kernel,npnts)).^2;

% initialize filtered signal matrix 
yFilt = zeros(1,length(Y));

% set initial conditions zi for filter function  
lead_dim = (max(length(1),length(filter_kernel))-1);
zi = zeros(1,lead_dim)';

% apply the filter to the data
[yFilt, zf] = filter(filter_kernel, 1, Y, zi, 2);
    
signalX3 = fft(yFilt);
signalAmp3 = 2*abs(signalX3)/npnts;

% plotting
xlimits = [20e3/1e3 120e3/1e3];
% freq response plot 
figure;
subplot(4,1,1)
stem(hz,signalAmp(1:length(hz)))
title('Frequency Response')
ylabel('Amplitude')
xlim(xlimits)

% freq response after filter
subplot(4,1,2)
stem(hz,signalAmp3(1:length(hz)))
ylabel('Amplitude')
xlim(xlimits)
title('Filtered Subband')

% ideal vs actual filter kernel
subplot(4,1,3)
plot(hz,filtkernX(1:length(hz)),'-','linew',2,'markersize',1)
hold on 
plot([0 passband(1) passband passband(2) nyquist]./1e3,[0 0 1 1 0 0],'ro-','linew',2,'markerfacecolor','w')
xlim(xlimits)
ylabel('Filter gain')
title('Frequency response of filter (fir1)')
legend({'Actual';'Ideal'}) 

%plot power
subplot(4,1,4)
plot(hz,10*log10(filtkernX(1:length(hz))))
xlim(xlimits)
xlabel('Frequency (Hz)')
ylabel('Filter gain (dB)')
title('Frequency response of filter (fir1)')

Code Output Plots:

Answer 1

From the comments the OP has further detailed the interest is in making a filter with a Gaussian frequency response. A Gaussian pulse in the time domain is also a Gaussian in the frequency domain. In general, the Fourier Transform of the impulse response is the frequency response, so if we seek a filter with a Gaussian frequency response, then it will have an impulse response also in the shape of a Gaussian. Extending this to discrete time: the frequency response is the Discrete Time Fourier Transform of the unit sample response, and for FIR filters, the unit sample response (discrete time impulse response) is the coefficients of the FIR filter.

To accomplish this I recommend using the windowing approach to filter design: determine the ideal impulse response (which here as a Gaussian pulse would have infinite support in time and is non-causal), sample that and then truncate and delay the sampled pulse for a causal implementation. To reduce the effects of truncation, multiply the resulting impulse response with a window. This will widen the resulting frequency response since the result in frequency will be the convolution of the Gaussian response with the window's Kernel, which can be compensated for by narrowing the target frequency response accordingly. If a bandpass is desired, the filter can then be translated to a real bandpass filter by modulating the coefficients with a sinusoid (demonstrated at the end).

The ideal Gaussian impulse response is given as:

$$g[n] = \frac{1}{\sigma\sqrt{2\pi}}e^{\frac{-n^2}{2\sigma^2}} \label{1}\tag{1}$$

Where $\sigma$ is the standard deviation in samples. We can related this to the 3 dB bandwidth of the frequency response using the Fourier Transform (The frequency response is the Discrete Time Fourier Transform, which given a long enough impulse response, this will reasonably approximate):

$$G(\omega) = e^{\frac{-\omega^2 \sigma^2}{2}}\label{2}\tag{2}$$

$G(\omega)$ is the normalized frequency response with frequency $\omega$ in units of radians\sample. Thus setting \ref{2} to $G(\omega_c)=1/\sqrt{2}$, we can determine $\sigma$ for a desired -3 dB bandwidth $\omega_c$:

$$\frac{1}{\sqrt{2}} = e^{\frac{-\omega_c^2 \sigma^2}{2}}$$

$$\sigma = \frac{\sqrt{ln(2)}}{\omega_c} \tag{3} \label{3}$$

Using this in equation \ref{1}, using a non-causal $n$ that extends symmetrically about $n=0$, we have the coefficients for the Gaussian filter with a 3 dB BW $\omega_c$ in radians/sample, prior to further windowing as needed to reduce sidelobes/ improve stopband rejection, as needed traded with the overall length of the filter.

Below is a demonstration for a 51 tap Gaussian FIR filter with a 3 dB cut-off of $\omega_c = 0.1$ radians/sample. With $\omega_c=0.1$, using Equation \ref{3} we get $\sigma = 8.32$ samples. With the computed $\sigma$, Equation \ref{1} is then used to compute the 51 filter coefficients with $n$ extending from -25 to +25.

The resulting coefficients and frequency response is plotted below. Note the floor in the frequency response due to the sharp truncation of the Gaussian pulse without further windowing. The orange curve in the frequency response is showing the true Gaussian response in frequency, as an estimate to what the actual frequency response of the Gaussian derived coefficients would be as given by the DTFT shown in blue. As the bandwidth is narrowed, the time response increases and the effects of such truncation become more pronounced, increasing the need for more coefficients and/or additional windowing.

By multiplying the derived coefficients with a Kaiser window (here with $\beta = 5$, the stopband is further pushed down and the bandwidth accordingly increased. If the enhanced stopband is required, my suggestion would be to simply iterate on the simulation as the target bandwidth is tightened to compensate (or the result can be derived mathematically from the convolution of the Kernels).

To convert the filter to a passband frequency at $\omega_1$, simply multiply the coefficients as computed with Equation \ref{1} with $2\cos(\omega_1 n)$. Below shows an example 151 tap bandpass design using a one-sided 3 dB BW $\omega_c = 0.03$ rad/sample at a band center $\omega_1=1 $ rad/sample.

---

Adding example MATLAB code at OP's request for a case with 250 KHz sampling rate, 10 KHz double-sideband 3dB bandwidth at a center frequency of 40 KHz. With 61 taps, no additional windowing is required, achieving rejection greater than 100 dB and exactly the 10 KHz BW. As bandwidth is narrowed, more taps will be needed to achieve similar performance (or windowing can be introduced in the total number of taps is an issue):

Zoom in confirming -3 dB at 35KHz to 45 KHz (10 KHz 3 dB BW):

MATLAB (Octave) Code below:

% Gaussian filter

fs = 250e3;  % sample rate in Hz
ntaps = 61;   % number of coeff (odd)
n = [-floor(ntaps/2): floor(ntaps/2)];

bw = 10e3;                   % desired 2 sided bw in Hz
fc = 40e3;                   % center frequency in Hz
wc = 2 * pi * bw/(2*fs);     % normalized one-sided bw rad/sample
sigma = sqrt(log(2))/wc;     % standard deviation for Gaussian in samples

# compute filter coefficients for lowpass Gaussian filter
coeff = 1/(sigma * sqrt(2*pi)) * exp(-n.^2/(2*sigma^2));

# convert lowpass filter to real bandpass
bpf = 2*coeff.*cos(2*pi*fc/fs*n);

# plot frequency response
[h,w] = freqz(bpf, 1, worN=2**14, fs=fs);
figure()
plot(w/1000, 20*log10(abs(h)))
grid()
axis([0, fs/2000, -150, 0])
xlabel("Frequency (KHz)")
ylabel("dB")

Adding additional frequency response for filter given by coeff below by replacing bpf with coeff in the freqz function above:

Answer 2

A solution I found myself is using the firpmord function in Matlab which selects the optimal filter order with bandwidth, amplitude min/max, error and sample rate inputs.

The order can then be plugged into the firpm function with the desired attenuation weighting to create a more optimized FIR filter.

I will say it is important to know what parameter you are trying to optimize for. Namely, attenuation inside and outside of filter or filter shape (avoiding rippling artifacts).

Here is the FIR kernel code that would need modified from my initial post:

% FIR filter specs
bandwidth = 20e3;
passband = [freq1-bandwidth freq1+bandwidth];
transw = 0.01;   %how the filter's edges taper

c = firpmord(passband,[90 0],[0.01 0.01],srate);  %optimize filter order for computation time

shape = [ 0 0 1 1 0 0 ];    %FIR shape

%define frequency shape of the FIR shape, firls function requires
%frequencies to go from 0 to 1
frex = [0 passband(1)-passband(1)*transw passband passband(2)+passband(2)*transw nyquist]/nyquist;

%define kernel
Wstop = 50;
Wpass = 1;

filter_kernel = firpm(c, frex, shape, [Wstop Wpass Wstop]);

% power spectrum of filter kernel
filtkernX = abs(fft(filter_kernel,npnts)).^2;

In my experience it seems the Gaussian is slightly better at creating a consistent amplitude for the filter over a wide variety of bandwidths while the firpm function requires strict amplitude settings in firmpord.