Digital Filter Design for Dummies
dmitry | Jan. 15, 2021, 6:39 p.m.
Author's Warning: I know little about digital filters except the absolute basics (and I may be wrong there as well).
A digital filter goal is to remove a signal's frequency components that outside's the filter's specifications (e.g. lowpass removes all frequency components above the filter's frequency limit, highpass removes all frequency components below the filter's frequency limit and bandpass removes all frequency components outside a frequency range). A filter itself can be perceived as an abstract series of numbers (my filters have the specific property of containing a symmetrical series of numbers) , where these numbers are known as coefficients. These coefficients are also known as taps (another term that is often used is the filter order, which is the number of taps - 1). This series of numbers is used in a mathematical operation with an input signal to produce an output signal without the specified frequency components. In my case, I'm looking to remove audible frequency components of a sound signal (where the signal is being caused by air pressure disturbances), leaving only the ultrasound components.
The signal itself can be anything, really not just sound signals. They key terms to remember are time domain and frequency domain. The data that most of us have usually worked with is in the time domain, i.e. think of stock price charts, which show stock price observations with time. A key theory I've been working with that is that any such signal in the time domain can be replaced by a series of sine curves with different amplitudes that are oscillating (alternating between their max and min amplitude values) at different frequencies (or number of times per second). Adding up this series of sine curves gives us the messy time domain charts we're used to working with. Visualizing that series of curves in the time domain by disassembling the overall signal in a chart will not be very useful to you. Instead, you can do a FFT on the time domain signal, which will give you back a chart showing the amplitudes and frequencies of the sine curves making up your signal. A filter's goal is to eliminate the sine curves with frequencies outside the filter specifications.
The effect of digital filter can be demonstrated by a chart showing how much the amplitude of sine curve at specific frequencies is reduced (meaning a response of 0 means no change). An amplitude decrease of -50 dB can represent that the signal at those frequencies is so small that it can be ignored.
In an ideal world, my bandpass filter looks like this:
The above chart shows that in a certain frequency range there is no change to amplitudes and that at the rest there is a reduction of -50 dB (dB is just a logarithmic system that represents a physical quantity with a very wide range, meaning that sensible single/double digit numbers represent both very small and very large numbers). The -50 dB is just an arbitrary way to represent a severe reduction in amplitude outside the passband of my bandpass (such a severe reduction is practically equivalent to eliminating those frequency sine components entirely).
However real bandpass filters have frequency response curves such as this:
The sharpness in the transition is completely lost and we have to deal with transition zones, where the frequency response curve slopes towards a target value of 0 or -50 (the passband area also loses flatness in its section of the frequency response curve as well, meaning our signal's amplitude at those frequencies will actually be slightly enhanced by different amounts, making the system's output less homogenous). The filters I'm trying to design also have an important design constraint; filter order. The filter order is the number of coefficients+1 that are used in the filtering mathematical operation. A higher filter order means that the resultant frequency response curve produced by the mathematical filtering operation has much more flexibility in how you want it to look, but conversely requires many more calculations when designing the filter or, more importantly, applying the filter to a signal.
This constraint comes from the specifications of hardware filters used in the project. These hardware filters can only store a set number of filter coefficients. If I was able to set a high filter order (~100) when designing the filter, the bandpass frequency response curve would look a lot closer to ideal (rather than the real band passes shown in the figure above).
Once you know the filter order, the type of bandpass and the frequency areas of interest (i.e. which frequencies you want to preserve/attenuate), you can design a filter, which really means use a fancy algorithm that takes your design inputs and outputs back filter design coefficients (i.e. the actual numbers used in the filtering operation). One additional design input you can play with is the width of the transition zones (i.e. how wide are those rolling curve sides on either end of the passband). When you're dealing with a filter order constraint, iterating these inputs becomes a game of tradeoff optimization (e.g. your passband becomes a little flatter but the transition zones becomes less steep and the overall attenuation in the stopbands decreases or there are more ripples that increase in size). I dealt exclusively with FIR filters that were symmetrical (this is the default filter type used in embedded microcontrollers).
Filter design can be done using Python's scipy library (the main FIR design functions they use are firls, firwin, firwin2 & remez. The freqz function is useful here as well since it calculates the frequency response of a given set of filter coefficients. The online documentation does a good job showing how to quickly make use of the functions, making me able to use some of Python's other modules to develop a little script that took a set of user inputs and displayed a chart of all the frequency responses (allowing you to quickly compare desired performance).
So first, I used PySimpleGUI to design the user input form (frequency values in kHz):
I then used mplcursors to plot the frequency responses of the 4 main FIR filter design algorithms to help compare what the filter performance would be:
Again, there is no perfect design algorithm.
Python is not the only tool you can use; when originally doing research on filter design, I found an excellent series of articles demonstrating filter design in Octave. While Octave is free and open source, MATLAB is my preferred choice due to the greater number of modules, functions, documentation etc. One excellent option available in MATLAB is an app called Filter Designer. It allows you to do what my Python script does, albeit in a cleaner GUI and with more filter design options (my Python script does allow you to compare multiple algorithms for the same inputs at the same time however!!). I would go so far to say that MATLAB's Filter Designer is perfect for a beginner trying to get an intuitive feel for the design tradeoffs inherent to filter design.
As an experiment to further investigate filter algorithm design tradeoffs, I created a MATLAB script that crunched through a variety of different input combinations based on the options that are most prominent in Filter Designer. This includes the passband range, stopband ranges, filter design algorithm and also two secondary inputs, density and weight (I wanted to confirm that these two inputs truly were secondary and had little influence on the overall frequency response curves at a low filter order constraint). It wasn't the most elegant script (not quite there yet) but by using a series of nested loops, I was able to vary all of the inputs while isolating changes to only one input at a time when comparing frequency responses. A set of inputs was then plugged in into a set of four different MATLAB filter design algorithms to calculate their frequency responses. Finally, rather than viewing thousands of charts (I did have to reduce the increment size for my input variables as the total number of combinations became quite large and the total computation time threatened to stretch into many days), I compared each calculated frequency response curve against an ideal bandpass and then saving any input selection that performed better than the previous best. I determined which performed "best" by measuring the root mean square error of each frequency response curve when compared with an ideal bandpass curve. Overall, the results of this experiment suggested that there is not much improvement that can be made with a fixed, low filter order (but remez seemed to perform best generally across a wide selection of desired bandpass ranges, except where this bandpass range becomes quite small, i.e. ~7kHz).
P.S Remez implementations are widely available online (Iowa Hills, Python scipy, MATLAB) but it is quite a complex algorithm.
P.P.S Actually "filtering" a signal or applying a filter's coefficients against an input is mathematically equivalent to convolution.
P.P.P.S Filtering can be done in several stages and in different ways (lowpass combined with highpass, two bandpasses together etc.). Order of filter stages does not matter.
Appendix: Code Samples
Python functions
#Function to design bandpass filters using Python's Remez, Firwin and firls algorithms
def des_bp_filt(des_inputs): #need to supply filter order, stop bands (low & high bounds), pass band (low & high bound) and sampling frequency
desired=(0,0,1,1,0,0)
order=int(des_inputs["Order"])
fs=int(des_inputs["Fs"])
bands=(int(des_inputs["SB1L"]), int(des_inputs["SB1H"]), int(des_inputs["PB1L"]), int(des_inputs["PB1H"]), int(des_inputs["SB2L"]), int(des_inputs["SB2H"]))
bands_firwin2=(0, int(des_inputs["SB1H"]), int(des_inputs["PB1L"]), int(des_inputs["PB1H"]), int(des_inputs["SB2L"]), fs/2)
bpass_coeffs_remez=signal.remez(order+1, bands,desired[::2], fs=fs)
bpass_coeffs_firls=signal.firls(order+1,bands, desired, fs=fs)
bpass_coeffs_firwin=signal.firwin(order+1, [int(des_inputs['PB1L']), int(des_inputs['PB1H'])],pass_zero=False, fs=fs)
bpass_coeffs_firwin2=signal.firwin2(order+1,bands_firwin2, desired, fs=fs)
return bpass_coeffs_remez, bpass_coeffs_firls, bpass_coeffs_firwin, bpass_coeffs_firwin2
#Function intending to allow for comparison of the four FIR design algorithms
def graph_filters(filt_dic, fs):
for key,coefficients in filt_dic.items():
freq, response = signal.freqz(coefficients, fs=int(fs))
plt.plot(freq,20* np.log10(abs(response)),label=key)
plt.legend(loc='best')
plt.title('Digital filter Frequency Responses\nClose chart window before selecting filter')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude (dB)')
plt.legend()
mplcursors.cursor()
plt.show()
#Primary window layout
layout_first = [[sg.Text("Step 1: Please fill in the below variables in Hz", size=(100, 2), justification='center', text_color='black')],[sg.Text("Stopband 1 Start (default is 0 "), sg.InputText(key='SB1L')],
[sg.Text("Stopband 1 End (default is 18000 "), sg.InputText(key='SB1H')], [sg.Text("Passband 1 Start (default is 20000 "), sg.InputText(key='PB1L')],
[sg.Text("Passband 1 End (default is 42000 "), sg.InputText(key='PB1H')],[sg.Text("Stopband 2 Start (default is 50000 "), sg.InputText(key='SB2L')],
[sg.Text("Stopband 2 End (default is Fs/2) "), sg.InputText(key='SB2H')], [sg.Text("Filter Order (default is 16) "), sg.InputText(key='Order')],
[sg.Text("Sampling Frequency (default is 253906)"), sg.InputText(key='Fs')], [sg.Button('Step 2: Design filters with chosen parameters. Close comparison chart before proceeding with Step 3', size=(100,2))],[sg.Text(" ", size=(100, 2))],
[sg.Text("Step 3: Please on your preferred filter below to start write process. ", size=(100, 2),text_color='black', justification='center')],
[sg.Listbox(filter_list, key='Select one filter only', size=(10, 10), enable_events=True)]]
#Confirmation coefficients were written to file
layout_second = [[sg.Text("Filter Coefficient write successful!")]]
window_first = sg.Window(title="Input Filter Design Variables and Select Preferred Filter", layout=layout_first)
while True:
event, window_values = window_first.read()
if event ==sg.WIN_CLOSED: # if user closes window or is satisfied with filters
window_first.close()
break
if event =='Step 2: Design filters with chosen parameters. Close comparison chart before proceeding with Step 3':
for key, value in window_values.items(): #Import values for non-empty fields
if value!="":
if key=="Fs":
design_params['SB2H']=int(value)/2
design_params[key]=value
remez, firls, firwin, firwin2 = des_bp_filt(design_params) #Calculate filter design coefficients
i=0
for fir in (remez, firls, firwin, firwin2):
filt_dict[filter_list[i]]=fir
print('{}\n{}'.format(filter_list[i],fir))
i=i+1
graph_filters(filt_dict, design_params['Fs']) #Graph all 4 designed filters
if event=='Select one filter only': #Select filter choice after reviewing and closing graph
filter_choice=window_values['Select one filter only'][0]
hex_coeffs=coeff_conv_hex(filt_dict[filter_choice])
write_coeffs_file(hex_coeffs)
trx_coeffs(hex_coeffs)
window_first.close()
window_second = sg.Window(title="Input Filter Design Variables and Select Preferred Filter", layout=layout_second, margins=(100, 50))
event, window_values = window_second.read()
if event==sg.WIN_CLOSED:
window_second.close()
MATLAB Script
function filt_opt(Fs, Order, Filt_LE, Filt_HE)
%initializing variables to store optimal filter settings
array_size=2048;
best_rho=0.5;
best_RMSE=100;
best_Fstop1=0;
best_Fstop2=0;
best_dens=0;
best_Wstop1=0;
best_Wstop2=0;
best_Wpass=0;
passband_width=Filt_HE-Filt_LE;
filter_index=0;
freq_xaxis=linspace(0,Fs/2,array_size);
bandpass_locations=find(freq_xaxis>Filt_LE & freq_xaxis<Filt_HE); %use logical indexing to find locations between the two bandpass ends
comparison_bandpass_locations=find(freq_xaxis>Filt_LE-(passband_width/6) & freq_xaxis<Filt_HE+(passband_width/6));
ideal_bandpass=zeros([1 array_size]); %bandpass area will aim for 0
ideal_bandpass(1:(bandpass_locations(1)-1))=-40; %Stopband 1 area will aim for -20
ideal_bandpass((bandpass_locations(end)+2):end)=-20; %Stopband 2 area will also aim for -20
ideal_bandpass=ideal_bandpass((comparison_bandpass_locations(1)):(comparison_bandpass_locations(end)));
for Fstop1=1000:1000:Filt_LE-1000
for Fstop2=Filt_HE+15000:-1000:Filt_HE+1000
%comparison_bandpass_locations=find(freq_xaxis<Fstop2); %using logical indexing to find locations below stopband 2. Assumptions is we care a lot about eliminating audible and having perfect flat response at 0 in the passband area. Don't care very much about performance in 40kHz+ area.
for dens=1:5:100
for Wstop1=0.2:0.2:1
for Wstop2=0.2:0.2:1
for Wpass=0.2:0.2:1
%tic
%firpm section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
firpm_bandpass_coeffs=firpm(Order, [0 Fstop1 Filt_LE Filt_HE Fstop2 Fs/2]/(Fs/2), [0 0 1 1 0 ...
0], [Wstop1 Wpass Wstop2], {dens});
firpm_lowpass_coeffs = firpm(Order, [0 Filt_HE Fstop2 Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop2], ...
{dens});
firpm_highpass_coeffs = firpm(Order, [0 Fstop1 Filt_LE Fs/2]/(Fs/2), [0 0 1 1], [Wstop1 Wpass], ...
{dens});
firpm_coeffs_x2=conv(firpm_bandpass_coeffs,firpm_bandpass_coeffs);
firpm_low_high_coeffs=conv(firpm_lowpass_coeffs,firpm_highpass_coeffs);
[h_1_1, ~]=freqz(firpm_bandpass_coeffs,1,array_size, Fs);
h_1_1=h_1_1(comparison_bandpass_locations(1):comparison_bandpass_locations(end));
[h_1_2, ~]=freqz(firpm_coeffs_x2,1,array_size, Fs);
h_1_2=h_1_2(comparison_bandpass_locations(1):comparison_bandpass_locations(end));
[h_1_3, ~]=freqz(firpm_low_high_coeffs,1,array_size, Fs);
h_1_3=h_1_3(comparison_bandpass_locations(1):comparison_bandpass_locations(end));
%firls section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
firls_bandpass_coeffs=firls(Order, [0 Fstop1 Filt_LE Filt_HE Fstop2 Fs/2]/(Fs/2), [0 0 1 1 0 ...
0], [Wstop1 Wpass Wstop2]);
firls_lowpass_coeffs = firls(Order, [0 Filt_HE Fstop2 Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop2]);
firls_highpass_coeffs = firls(Order, [0 Fstop1 Filt_LE Fs/2]/(Fs/2), [0 0 1 1], [Wstop1 Wpass]);
firls_coeffs_x2=conv(firls_bandpass_coeffs,firls_bandpass_coeffs);
firls_low_high_coeffs=conv(firls_lowpass_coeffs,firls_highpass_coeffs);
[h_2_1, ~]=freqz(firls_bandpass_coeffs,1,array_size, Fs);
h_2_1=h_2_1(comparison_bandpass_locations(1):comparison_bandpass_locations(end));
[h_2_2, ~]=freqz(firls_coeffs_x2,1,array_size, Fs);
h_2_2=h_2_2(comparison_bandpass_locations(1):comparison_bandpass_locations(end));
[h_2_3, ~]=freqz(firls_low_high_coeffs,1,array_size, Fs);
h_2_3=h_2_3(comparison_bandpass_locations(1):comparison_bandpass_locations(end));
%firgr section %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
firgr_bandpass_coeffs=firgr(Order, [0 Fstop1 Filt_LE Filt_HE Fstop2 Fs/2]/(Fs/2), [0 0 1 1 0 ...
0], [Wstop1 Wpass Wstop2], {dens});
firgr_lowpass_coeffs = firgr(Order, [0 Filt_HE Fstop2 Fs/2]/(Fs/2), [1 1 0 0], [Wpass Wstop2], ...
{dens});
firgr_highpass_coeffs = firgr(Order, [0 Fstop1 Filt_LE Fs/2]/(Fs/2), [0 0 1 1], [Wstop1 Wpass], ...
{dens});
firgr_coeffs_x2=conv(firgr_bandpass_coeffs,firgr_bandpass_coeffs);
firgr_low_high_coeffs=conv(firgr_lowpass_coeffs,firgr_highpass_coeffs);
[h_3_1, ~]=freqz(firgr_bandpass_coeffs,1,array_size, Fs);
h_3_1=h_3_1(comparison_bandpass_locations(1):comparison_bandpass_locations(end));
[h_3_2, ~]=freqz(firgr_coeffs_x2,1,array_size, Fs);
h_3_2=h_3_2(comparison_bandpass_locations(1):comparison_bandpass_locations(end));
[h_3_3, ~]=freqz(firgr_low_high_coeffs,1,array_size, Fs);
h_3_3=h_3_3(comparison_bandpass_locations(1):comparison_bandpass_locations(end));
h_array=[20*log10(abs(h_1_1)), 20*log10(abs(h_1_2)), 20*log10(abs(h_1_3)), 20*log10(abs(h_2_1)),...
20*log10(abs(h_2_2)), 20*log10(abs(h_2_3)), 20*log10(abs(h_3_1)),...
20*log10(abs(h_3_2)), 20*log10(abs(h_3_3))]; %place all frequency response arrays into one array
%summary_inputs=[Fstop1; Fstop2; dens; Wstop1; Wstop2; Wpass];
%disp(summary_inputs);
for fir_filter=1:1:size(h_array,2)
y_yhat_sq2=(transpose(h_array(:,fir_filter))-ideal_bandpass).^2;
RMSE = sqrt(mean(y_yhat_sq2, 'all'));
pearson_coeff_rho=corrcoef(ideal_bandpass,transpose(h_array(:,fir_filter)));
if abs(RMSE)<abs(best_RMSE)
best_RMSE=RMSE;
best_rho=pearson_coeff_rho(2);
best_Fstop1=Fstop1;
best_Fstop2=Fstop2;
best_dens=dens;
best_Wstop1=Wstop1;
best_Wstop2=Wstop2;
best_Wpass=Wpass;
%best_filt_coeffs=coeff_array(t,:);
filter_index=fir_filter;
%disp(best_rho);
end
end
%toc
end
end
end
end
end
end
fileID = fopen('Optimal_filters.txt','a+');
fprintf(fileID,'For inputs of Filt_LE %i, Filt_HE %i, Order %i and Fs %i\n', Filt_LE, Filt_HE, Order, Fs);
fprintf(fileID,'Optimal filter is selection %i, best_Fstop1 is %i, best_Fstop2 is %i, best_dens is %i, best_Wstop1 is %i, best_Wstop2 is %i, best_Wpass is %i. This has a rho of %i and a RMS of %i\n', ...
filter_index, best_Fstop1, best_Fstop2, best_dens, best_Wstop1, best_Wstop2, best_Wpass, best_rho, best_RMSE);
%fprintf('Filter Coefficients: %i', best_filt_coeffs);
fclose(fileID);
end
About Us
I find answers for questions nobody else has time to answer and to help me remember, I write them dowm here!
0 comments