Optimize FIR filter kernel in Matlab

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

% 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:

• I am confused by some of your statements— if the ideal shape you want in time is Gaussian, then the Kernel in frequency will also be Gaussian. Your ideal Kernel that you show in frequency as a brickwall filter would be a Sinc in time. Am I missing something in my interpretation of what you are doing or what you want to do? Commented Apr 10, 2023 at 19:38
• @DanBoschen I want a gaussian-looking filter kernel in frequency domain as shown in the left-most plot above and the third subplot. I want a smooth filter that isn't distorted at its peak. With my current integration approach I am getting distortion at the peak as shown on the right-most plot. I want a Gaussian-like filter with a distinct peak at a gain of 1. My question is how do I create an algorithm that picks the order for this kernel that gives theses characteristics given a range of frequencies I want to filter. Commented Apr 10, 2023 at 20:28

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:

• I'm confused with what you are showing. I already have the filter defined in matlab via the firls function. I just need a way to make sure the kernel doesn't get distorted when I am selecting the order Commented Apr 11, 2023 at 16:45
• @PrematureCorn your issue is using firls to create a Gaussian filter. I am showing you an alternate approach that will result in a much better match to the filter shape you want. Commented Apr 11, 2023 at 16:54
• (Given you know you want a Gaussian response, and we know what that function is, you needn’t compute a least squares estimate to such response- we already have the result directly from the Gaussian function— try it, I think it will do what you want very nicely) Commented Apr 11, 2023 at 17:28
• I am happy to try your approach, I am a bit confused on the kaiser multiplication and truncated gaussian. I was expecting a bandpass response in dB scale similar to what I am showing in my last subplot. Could you share some of your matlab so I could play around with your approach and understand better? Commented Apr 11, 2023 at 18:06
• I'd like to accommodate bandwidths between 30kHz and 2kHz. A nice test center freq would be 40kHz and attenuation of -30dB Commented Apr 11, 2023 at 18:39

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.

• firpm is the Parks McLellan algorithm for creating optimum filters in the min max sense (minimizes peak error) rather than least squares error such as firls. In most applications I use, having a minimum least squares error is superior (try yourself with same order and compare firls to firpm-- firls will have much less error overall (in an rms or least squares sense) while firpm will have less peak error - which occurs right at the band edge. In any event, my understanding was you wanted a Gaussian shaped frequency response (as I showed in my answer), or was that not really your objective? Commented Apr 15, 2023 at 2:45