Assistance Needed with FIR Comb Peak and Notch Filter Design in MATLAB

Hello,

I am working with a sampling frequency Fs=10Fs = 10Fs=10 MHz and a fundamental frequency of 0.5 MHz. I need to design two types of FIR comb filters in MATLAB, but I’m unsure how to write the code for these filters manually. Previously, I used IIR filters generated with filterDesigner, but I encountered phase distortion issues, which may lead to errors in my data analysis. I now want to switch to FIR filters to avoid this phase distortion.

Here’s what I’m looking to implement:

1. Comb Peak Filter (1-3 MHz): This filter should capture higher harmonic components at 1 MHz, 1.5 MHz, 2 MHz, 2.5 MHz, and 3 MHz, each with a bandwidth of 0.1 MHz.
2. Comb Notch Filter (3-4 MHz): This filter should isolate four wideband noise bands (not including subharmonic or superharmonic components) centered around:

• 3-3.2 MHz
• 3.3-3.45 MHz
• 3.55-3.7 MHz
• 3.8-4 MHz

I’d like to create these FIR filters but am unsure how to set up the design parameters or write the code manually. Any guidance or sample code on how to implement these filters as FIR filters would be greatly appreciated.

Thank you very much for your help!

[Mathieu NOE]

hello

you want to switch to FIR filters to have a linear phase vs frequency ?

[轶凡]

Hello, yes, I am looking to switch to FIR filters primarily to achieve a linear phase response across frequencies. My current IIR filters introduce phase distortion, which affects my subsequent data processing. But I am currently using the filtfilt function, and I am not sure if this can eliminate the phase distortion of the IIR filter.

[William Rose]

@轶凡, yes filtfilt() will eliminate the phase distroption of an IIR filter.

If you want to use IIR filters and filtfilt, then, for the comb peak filter, you can make 4 bandpass filters and pass the input through all 4 in pranllel, then add the outputs to get the combined output:

y=filt1(x)+filt2(x)+filt3(x)+filt4(x)

where filt1(), etc. represents an IIR filter implemented with filtfilt().

The trouble I anticipate with this plan is that you may need a high order IIR filter to acheive the narrow bands you desire. filtfilt() requires you to work with the [b,a] specification of filters. Unfortunately the [b,a] representation of filters is numerically less stable than second order secotions or z,p,k. So you may not be able to get the performance you want with filtfilt(). If so, then you will use FIR. You will have to experiment to find out.

[轶凡]

Thank you for your answer. Below, I have shown the amplitude and phase response results of my four IIR and FIR filters. From the graphs, it seems that the phase response of the FIR filters fluctuates more. As a beginner in filtering, I learned that using FIR filters is one way to solve phase distortion. I may not fully understand the phenomenon of fluctuating FIR phase response.

Answer 1

[Star Strider]

The peak filter was a bit challenging, however I got it to work reasonably well.

The comb filter may not exactly meet your design specifications because the stopbands are very narrow and very close together.

This is the best I can do with your filters —

Fs = 1E+7;
fv = (1:0.5:3)*1E+6;
fm = fv(:) + [-1.25 -1 0 1 1.25]*1E+5 + 5E+4;                           % Filter Stopbands
BW = fm(:,4) - fm(:,3)                                                  % Design Bandwidths

BW = 5×1

100000 100000 100000 100000 100000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fcomb = reshape(fm.', 1, []);
fcomb(end) = [];
% mags = [[1 0 1] reshape((zeros(numel(fcomb)/4-2,1)+[0 1]).', 1, [])  [0 1]];
% dev = [[0.5 0.1 0.5] reshape((zeros(numel(fcomb)/4-2,1)+[0.1 0.5]).', 1, []) [0.1 0.5]];
mags = [[0 1 0] reshape((zeros(numel(fcomb)/4-2,1)+[1 0]).', 1, [])  [1 0]];
dev = [[0.1 0.5 0.1] reshape((zeros(numel(fcomb)/4-2,1)+[0.5 0.1]).', 1, []) [0.5 0.1]];
[n,Wn,beta,ftype] = kaiserord(fcomb,mags,dev,Fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');

fprintf('Filter Order = %4d so the signal must have at least %6d elements', n, 3*n)

Filter Order =  336 so the signal must have at least   1008 elements

figure
freqz(hh, 1, 2^22, Fs)
sgtitle('Peak Filter')

BW = [3 3.2; 3.3 3.45; 3.55 3.7; 3.8 4]*1E6                             % Comb Filter Stopbands

BW = 4×2

3000000 3200000 3300000 3450000 3550000 3700000 3800000 4000000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fm = [0.99*BW(:,1) BW 1.01*BW(:,2)]

fm = 4×4

2970000 3000000 3200000 3232000 3267000 3300000 3450000 3484500 3514500 3550000 3700000 3737000 3762000 3800000 4000000 4040000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

BW = fm(:,4) - fm(:,3)                                                  % Design Bandwidths

BW = 4×1

32000 34500 37000 40000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% fv = (1:0.5:3)*1E+6
% fm = fv(:)+[-5 -2 2 5];
fm = reshape(fm.',1,[]);
% fcomb = reshape(fm.', 1, []);
% fcomb(end) = [];
fcomb = fm

fcomb = 1×16

2970000 3000000 3200000 3232000 3267000 3300000 3450000 3484500 3514500 3550000 3700000 3737000 3762000 3800000 4000000 4040000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

mags = [[1 0 1] reshape((zeros(numel(fcomb)/4-2,1)+[0 1]).', 1, [])  [0 1]];
dev = [[0.5 0.1 0.5] reshape((zeros(numel(fcomb)/4-2,1)+[0.1 0.5]).', 1, []) [0.1 0.5]];
[n,Wn,beta,ftype] = kaiserord(fcomb,mags,dev,Fs);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');

fprintf('Filter Order = %4d so the signal must have at least %6d elements', n, 3*n)

Filter Order =  280 so the signal must have at least    840 elements

figure
freqz(hh, 1, 2^22, Fs)
sgtitle('Comb Filter')

A better approach to the comb filter could be running four elliptic IIR filters in parallel.

.

[轶凡]

Thank you for your response！ I used four elliptical IIR bandpass filters to implement the first functionality. Below is the source code, along with the amplitude and phase responses of the four IIR filters. I'm not entirely sure if filtfilt() can achieve zero phase.

Fs = 10e6;
center_freqs = [1e6, 1.5e6, 2e6, 2.5e6];
bandwidth = 0.1e6;
for i = 1:length(center_freqs)
    f_low = (center_freqs(i) - bandwidth/2) / (Fs/2);
    f_high = (center_freqs(i) + bandwidth/2) / (Fs/2);
    [n, Wn] = ellipord([f_low, f_high], [f_low - 0.05, f_high + 0.05], 1, 60);
    [b, a] = ellip(n, 1, 20, Wn, 'bandpass');
    figure(i);
    [H, f] = freqz(b, a, Fs);
    subplot(2, 1, 1);
    hold on;
    plot(f / 1e6, 20 * log10(abs(H))); 
    title('single filter amplitude response');
    xlabel('frequency (MHz)');
    ylabel('amplitude (dB)');
    grid on;subplot(2, 1, 2);
    hold on;
    plot(f / 1e6, angle(H) * 180 / pi); 
    title('single filter phase response');
    xlabel('frequency (MHz)');
    ylabel('phase (degree)');
    grid on;
    % y_filtered = filtfilt(b, a, x);
    % y = y + y_filtered;
end

%At the same time, I also tried using FIR filters, but it seems that FIR might not perform the filtering task effectively.
Fs = 10e6;
center_freqs = [1e6, 1.5e6, 2e6, 2.5e6];
bandwidth = 0.1e6;
% y = zeros(size(rf_matrix,1),1);
for i = 1:length(center_freqs)
    f_low = (center_freqs(i) - bandwidth/2) / (Fs/2);
    f_high = (center_freqs(i) + bandwidth/2) / (Fs/2);
    order = 200;
    b = fir1(order, [f_low,f_high],'bandpass',kaiser(order+1,8));
    delay = round(order / 2);
    % y_filtered = filter(b, 1, mean(rf_matrix, 2));
    % y_filtered = [y_filtered(delay+1:end); zeros(delay, 1)];
    % y = y + y_filtered;
    % [step, Wn] = ellipord([f_low, f_high], [f_low - 0.05, f_high + 0.05], 1, 60);
    % [b, a] = ellip(step, 1, 60, Wn, 'bandpass');
    figure(i+4);
    [H, f] = freqz(b, 1, 2^14, Fs);
    subplot(2, 1, 1);
    hold on;
    plot(f / 1e6, 20 * log10(abs(H))); 
    title('single filter amplitude response');
    xlabel('frequency (MHz)');
    ylabel('amplitude (dB)');
    grid on;    
    subplot(2, 1, 2);
    hold on;
    plot(f / 1e6, angle(H) * 180 / pi); 
    title('single filter phase response');
    xlabel('frequency (MHz)');
    ylabel('phase (degree)');
    grid on;    
    % y = y + y_filtered;
end

[Star Strider]

My pleasure!

I though my FIR filters did reasonably well! In any event, they did about as well as FIR fiilters can do in such situations.

The filtfilt function will come very close to achieving zero phase. (I have actually never seen it fail.) Its only problem can be initial transients with a signal that does not begin with a zero value, and they can be dealt with by pre-appending a series of values (the length being about 100 times the sampling interval) that are the same amplitude as the first value of the signal, filtering that, and then removing that initial vector from the filtered signal. This is not a common problem with filtfilt, however if it is one with a particular signal, that is the way I usually deal with it.

A better way to iimplement your elliptic filters would be to change this:

[b, a] = ellip(n, 1, 20, Wn, 'bandpass');

to:

[z, p, k] = ellip(n, 1, 20, Wn, 'bandpass');

[sos,g] = zp2sos(z, p, k);

and then use those outputs with filtfilt to filter your signal. The second-order-section realisation generally produces a more stable filter and a more reliable fiiltered result. Alternatively, you could skip the filter design steps altogether and use the bandpass function, with the 'ImpulseResponse','iir' name-value paiir to get a reliable elliptic second-order-section filter realisation. Get both outputs, and use the second with filtfilt in subsequeent filterinig operations so that you don’t have to re-design the same filter each time.

If you want to create single-passband fiilters, it is generally preferable to use IIR filters (I prefer elliptic filters for their computational efficiency) rather than FIR filters, since FIR fiilters are best (most efficiient) for multiple passband or stopband filters, and may fail with signals that are too short in length. Use elliptic filters instead in those applicatiions.

.

[轶凡]

Thank you again for your response and suggestions. Here is the function and filtering results that I modified based on your advice.

function fil_rf_matrix = Notch_10MHz(rf_matrix, Fs)
    % this function aim to generate the comb filters using ellipse IIR filter
    % input:rf_matrix(1664*1024),samples*channels 
    % Fs:sample frequency
    % output:fil_rf_matrix(1664*1024)
    
    pre_length = 100*10;
    pre_values = repmat(rf_matrix(1,:), pre_length, 1);
    rf_matrix_line_exp = [pre_values; rf_matrix];
    fil_rf_matrix = zeros(size(rf_matrix));
    center_freqs = [1e6, 1.5e6, 2e6, 2.5e6];
    bandwidth = 0.1e6;
    for no_fil = 1:length(center_freqs)
        f_low = (center_freqs(no_fil) - bandwidth/2) / (Fs/2);
        f_high = (center_freqs(no_fil) + bandwidth/2) / (Fs/2);
        [step, Wn] = ellipord([f_low, f_high], [f_low - 0.05, f_high + 0.05], 1, 60);
        [z,p,k] = ellip(step, 1, 60, Wn, 'bandpass');
        [sos, g] = zp2sos(z,p,k);
        filtered_matrix = filtfilt(sos, g, rf_matrix_line_exp);
        fil_rf_matrix_temp = filtered_matrix(pre_length+1:end,:);
        fil_rf_matrix = fil_rf_matrix+fil_rf_matrix_temp;
    end
end

[Star Strider]

As always, my pleasure!

I have no ideea what you are filtering, however I am very happy that everything worked!