Home > matlab code, Uncategorized, youtube demo code > Filter Design Using Matlab Demo

## Filter Design Using Matlab Demo

This is the code I use in my youtube videos on filtering (ddorran channel)

%generate the noisy signal which will be filtered
x= cos(2*pi*12*[0:0.001:1.23]);
x(end) = [];
[b a] = butter(2,[0.6 0.7],'bandpass');
filtered_noise = filter(b,a,randn(1, length(x)*2));
x = (x + 0.5*filtered_noise(500:500+length(x)-1))/length(x)*2;
% This script is available at https://dadorran.wordpress.com search for
% filtering matlab demo
plot(x)
title('Noisy signal')
xlabel('Samples');
ylabel('Amplitude')

%plot magnitude spectrum of the signal
X_mags = abs(fft(x));
figure(10)
plot(X_mags)
xlabel('DFT Bins')
ylabel('Magnitude')

%plot first half of DFT (normalised frequency)
num_bins = length(X_mags);
plot([0:1/(num_bins/2 -1):1], X_mags(1:num_bins/2))
ylabel('Magnitude')

%design a second order filter using a butterworth design technique
[b a] = butter(2, 0.3, 'low')

%Show the signal flow diagram of the filter/system
create_signal_flow(b,a)

%plot the frequency response (normalised frequency)
H = freqz(b,a, floor(num_bins/2));
hold on
plot([0:1/(num_bins/2 -1):1], abs(H),'r');

%filter the signal using the b and a coefficients obtained from
%the butter filter design function
x_filtered = filter(b,a,x);

%plot the filtered signal
figure(2)
plot(x_filtered,'r')
title('Filtered Signal - Using Second Order Butterworth')
xlabel('Samples');
ylabel('Amplitude')

%Redesign the filter using a higher order filter
[b2 a2] = butter(20, 0.3, 'low')

%Plot the magnitude spectrum and compare with lower order filter
H2 = freqz(b2,a2, floor(num_bins/2));
figure(10)
hold on
plot([0:1/(num_bins/2 -1):1], abs(H2),'g');

%filter the noisy signal and plot the result
x_filtered2 = filter(b2,a2,x);
figure(3)
plot(x_filtered2,'g')
title('Filtered Signal - Using 20th Order Butterworth')
xlabel('Samples');
ylabel('Amplitude')

%Use a band reject filter in place of a low pass filter
[b_stop a_stop] = butter(20, [0.5 0.8], 'stop');

%plot the magnitude spectrum
H_stopband = freqz(b_stop,a_stop, floor(num_bins/2));
figure(10)
hold on
plot([0:1/(num_bins/2 -1):1], abs(H_stopband),'c');

%plot filtered signal
x_filtered_stop = filter(b_stop,a_stop,x);
figure(4);
plot(x_filtered_stop,'c')
title('Filtered Signal - Using Stopband')
xlabel('Samples');
ylabel('Amplitude')

%Use matlabs built-in buttord function to get the optimum order to meet a specification
[N Wn] = buttord(0.1, 0.5, 5, 40)

%use the N and Wn values obtained above to design the filter in the usual way
[b3 a3] = butter(N, Wn, 'low');

%plot the magnitude spectrum
H3 = freqz(b3,a3, floor(num_bins/2));
figure(10);
hold on
plot([0:1/(num_bins/2 -1):1], abs(H2),'k');
figure(10)

%filter the signal and plot the ouput of the filter
x_filtered3 = filter(b3,a3,x);
figure(5);
plot(x_filtered3,'k')
title(['Filtered Signal - Using ' num2str(N) ' th Order Butterworth'])
xlabel('Samples');
ylabel('Amplitude')
% comparison with other filter design techniques (chebyshev and elliptical)
[b_butter a_butter] = butter(4, 0.2, 'low');
H_butter = freqz(b_butter, a_butter);

[b_cheby a_cheby] = cheby1(4, 0.5, 0.2, 'low');
H_cheby = freqz(b_cheby, a_cheby);

[b_ellip a_ellip] = ellip(4, 0.5, 40, 0.2, 'low');
H_ellip = freqz(b_ellip, a_ellip);

%plot each filter to compare
figure(11)
norm_freq_axis = [0:1/(512 -1):1];
plot(norm_freq_axis, abs(H_butter))
hold on
plot(norm_freq_axis, abs(H_cheby),'r')
plot(norm_freq_axis, abs(H_ellip),'g')
legend('Butterworth', 'Chebyshev', 'Elliptical')
xlabel('Normalised Frequency');
ylabel('Magnitude')

%plot in dB for verification that spec is met
figure(12);
plot(norm_freq_axis, 20*log10(abs(H_butter)))
hold on
plot(norm_freq_axis, 20*log10(abs(H_cheby)),'r')
plot(norm_freq_axis, 20*log10(abs(H_ellip)),'g')
legend('Butterworth', 'Chebyshev', 'Elliptical')
xlabel('Normalised Frequency ');
ylabel('Magnitude (dB)')


1. November 20, 2014 at 9:26 pm

Thanks dude. It’s a great help.

2. January 11, 2015 at 4:42 pm

it is possible to create a simple low pass filter like in RC circuits? For instance if we create a sine wave like y=10*sin(2*pi*f*t). Can i just filter the signal with a cutoff frequency to be able to see how the filtered signal looks like(for the amplitude decay)?
so when i plug in some information for example f= cutoff and plot the filtered signal this must be reduced in amplitude somewhere at 30%percent of the input signal ?

3. January 11, 2015 at 7:46 pm

Thank you!
I am using butter filter design for noise removal from a voice signal and your videos have been of great help. Is there a way to design coefficients for the filter function which include more than a stop band? I would like to ‘compose’ a complicated impulse response and then only do one filtering operation, instead of filtering-designing-filtering…

• January 13, 2015 at 6:18 am

You could design your system using frequency response using pole-zero placement (see https://www.youtube.com/watch?v=m5TP2uG_O2M). You can then determine the transfer function H(z) from the positions of the poles and zeros.

4. February 24, 2015 at 7:31 pm

Sir, your videos and codes are of great help.But could you please explain the mapping of DFT Bin scale to Normalized frequency scale that you used in your code??

• February 24, 2015 at 7:44 pm

Normalised frequency maps a value of 1 to Nyquist Frequency (fs/2) Hz i.e. half the sampling frequency fs Hz.
A normalised frequency value of 0.5 therefore maps to fs/4 Hz;
A normalised frequency value of 0 maps to 0 Hz.

Bin number N/2 also corresponds to Nyquist frequency, where N is the total number of bins.
A normalised frequency of 1 therefore maps to bin number N/2.
A normalised frequency value of 0.5 maps to bin number N/4.
A normalised frequency value of 0 maps to bin number 0.

5. July 26, 2016 at 10:54 pm