Designing an Arbitrary FIR Filter

I need to design an Arbitrary FIR Filter whose response looks like this,

Frequency = [100 125 160 200 250 315 400 500 630 800 1000 1250 1600 2000 2500 3150 4000 5000 6300 8000 10000]

Required_Filter_Response = [55,38 55,74 52,13 52,64 51,25 49,04 46,50 44,99 42,95 40,73 39,47 37,53 35,36 34,33 32,33 30,36 28,01 25,98 23,19 20,42 18,65]

My professor suggested me using designfilt and fvtool for realizing it. I am not new to MATLAB, but to filters and digital signal processing, I am.

Any suggestions on how to achieve it would be greatly appreciated. Thanks!

Answer 1

[Star Strider]

This is probably as close as you can get:

Frequency = [100  125  160  200  250  315  400  500  630  800  1000  1250  1600  2000  2500  3150  4000  5000  6300  8000  10000];
Required_Filter_Response = [55.38  55.74  52.13  52.64  51.25  49.04  46.50  44.99  42.95  40.73  39.47  37.53  35.36  34.33  32.33  30.36  28.01  25.98  23.19  20.42 18.65];
Fs = 20000;                                                                 % Create Sampling Frequency (Not Supplied, Must Be > 20000 Here)
Fn = Fs/2;
Fnew = linspace(min(Frequency), max(Frequency), fix(numel(Frequency)/2)*2);
Rnew = interp1(Frequency, Required_Filter_Response, Fnew);
b = firls(100, Fnew/Fn, Rnew/max(Rnew));

figure
freqz(b, 1, 2^14, Fs)

[h,f] = freqz(b, 1, 2^14, Fs);
figure
plot(Fnew, 20*log10(Rnew/max(Rnew)))
hold on
plot(f, 20*log10(abs(h)))
hold off

Note that the filter and the desired response are essentially the same in the second plot.

[Krishna Ganesan]

Thank you very much for the quick response and the code.

My sampling frequency comes from reading in an audio file which is 44100 Hz. I have to apply this filter on the audio file.

However, when I run the script changing the sampling frequency (Fs) from 20000 Hz to 44100 Hz, MATLAB gives a warning.

and the plots don't look like the desired response.

Kindly help me find the issue.

[Star Strider]

My pleasure.

It is always good to know the sampling frequency annd the intended application.

That problem is then relatively straightforward to solve. It helps to estimate the parameters of mathematical model of the frrequency response, in order to see what the correct extrapolation values would be, out to the Nyquist frequency. This appears to be a reasonable model:

fcn = @(b,f) b(1).*f.^b(2) + b(3);
B = fminsearch(@(b) norm(Required_Filter_Response - fcn(b,Frequency)), [25; -0.15; 10])

producing

B =
          345.886358298793
       -0.0285538321269912
         -245.322882425467

So evaluating the function at the Nyquist frequency (22050 Hz) and adding those to the original vectors:

Frequency = [100  125  160  200  250  315  400  500  630  800  1000  1250  1600  2000  2500  3150  4000  5000  6300  8000  10000 22050];
Required_Filter_Response = [55.38  55.74  52.13  52.64  51.25  49.04  46.50  44.99  42.95  40.73  39.47  37.53  35.36  34.33  32.33  30.36  28.01  25.98  23.19  20.42 18.65 14.57];
Fs = 44100;                                                                 % Sampling Frequency
Fn = Fs/2;
Fnew = linspace(min(Frequency), max(Frequency), fix(numel(Frequency)/2)*2);
Rnew = interp1(Frequency, Required_Filter_Response, Fnew);
b = firls(60, Fnew/Fn, Rnew/max(Rnew));

figure
freqz(b, 1, 2^14, Fs)

[h,f] = freqz(b, 1, 2^14, Fs);
figure
plot(Fnew, 20*log10(Rnew/max(Rnew)))
hold on
plot(f, 20*log10(abs(h)))
hold off

So with the added frequency point and a slightly different (and more efficient) filter order, the new transfer function appears to almost exactly match ‘Required_Filter_Response’. It appears to be a linear-phase filter, however I would still recommend using it with filtfilt rather than filter.

EDIT — (17 Feb 2020 at 17:25)

Added Bode plot figure —

and code creating it —

[h,f] = freqz(b, 1, 2^14, Fs);
figure
subplot(2,1,1)
plot(Fnew, 20*log10(Rnew/max(Rnew)),'.r')
hold on
plot(f, 20*log10(abs(h)),'-b')
hold off
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
xlim([min(f) max(f)])
legend('Desired Filter Response','Actual Filter Response')
grid
subplot(2,1,2)
plot(f, rad2deg(unwrap(angle(h))))
xlabel('Frequency (Hz)')
ylabel('Phase (degrees)')
xlim([min(f) max(f)])
grid

This is simply an augmented freqz plot. It is not part of the filter design code.

[Krishna Ganesan]

Dear Star Strider,

Thank you very much for providing the code for the sampling frequency that I wanted. Still, after running and studying the code, I have the following doubts.

--------------------------------------------------------------------------------------------------------------------------------------------------

Doubt 1 - The filter response that I need is from 55 dB at 100 Hz to 20 dB at 10,000 Hz. I understand that the frequency is normalized using,

plot(Fnew, 20*log10(Rnew/max(Rnew)))

and I understand that normalized frequencies are always adjusted from 0 to 1. But, when I run the code, the filter response is from 0 dB at 0 Hz to -12 dB at 20,000 Hz as presented below,

What my doubt is, should some sort of scaling be done? Or any gain be added to boost the filter response to make it to the scaling (from 55 dB at 100 Hz to 20 dB at 10,000 Hz) that I want?

--------------------------------------------------------------------------------------------------------------------------------------------------

Doubt 2 - I have an audio file with two channels on which the filter has to be applied. I read in the audio file using,

[y,Fs]    = audioread('Input.wav'); %Audio measurement recorded using Binaural headset microphones
Left_Ear  = y(:,1);
Right_Ear = y(:,2);

My doubt here is, how do I apply the filter on the audio file?

I understand that the filter is stored in the variable 'b'

b = firls(60, Fnew/Fn, Rnew/max(Rnew));

Would this be a correct approach to apply the filter on the audio file? As my professor wanted me to apply the filter in time-domain.

Filtered_Left_Ear  = conv(Left_Ear,b);
Filtered_Right_Ear = conv(Right_Ear,b);

I understand that the filter has to be applied to each channels separately and then they have to be appended together to save the filtered audio file.

Filtered_Audio_File = [Filtered_Left_Ear, Filtered_Right_Ear];
audiowrite('Output.wav', Filtered_Audio_File, Fs);

--------------------------------------------------------------------------------------------------------------------------------------------------

Kindly help me solve the issue

[Star Strider]

My pleasure.

The filter is standardised such that the maximum amplitude is considered to be 0 dB, with the frequency defined between 0 Hz and the Nyquist frequency. That is the convention. The filter works as designed.

If you want to amplify the input or output by a scalar constant, the just multiply it. The filter will continue to work as designed.

Otherwise, if you want to re-design the filter to other specifications, you are certainly welcome to do so. I designed it to the specifications you posted.

Filter the signal in the time domain as:

Filtered_Left_Ear = filtfilt(b, 1, Left_Ear);
Filtered_Right_Ear = filtfilt(b, 1, Right_Ear);

Since ‘Left_Ear’ and ‘Right_Ear’ are identically-size columns in the same matrix, you can use that matrix as input to filtfilt, that will then return a matrix of the same size:

y = [Left_Ear(:), Right_Ear(:)];
Filtered_Ears = filtfilt(b, 1, y);

That is all there is to it.

[Krishna Ganesan]

Thank you very much for the clarification!

[Star Strider]

My pleasure.

If my Answer helped you solve your problem, please Accept it!