FIR Low-Pass Filter Design with Window Method

Load Signal and Parameters

We load the signal we need.

load('tecla_piano.mat');

We calculate the sampling period and the number of samples.

Ts = 1/fs;
N = length(data);
n = 0:N-1;
nTs = n * Ts;

Set the number of frequency points which generate the spectrums. We also generate the frequency vector.

nFFT = 2 * fs;
f = (0:nFFT-1) * fs/(nFFT-1);

Time Domain Signal Representation

We represent the input signal in the time domain.

figure;
plot(nTs,data);
title('Input signal, TIME domain');
xlabel('nT_{S} (s)');
ylabel('data[nT_{S}] (V)');
grid on;

Signal Spectrum Analysis

We calculate and observe the spectrum of the signal. From the magnitude spectrum, we decide that we want to filter the signal up to 1000 Hz and we want to remove it completely from 1500 Hz up.

data_f = fft(data,nFFT) / min(N,nFFT);
figure;
subplot(2,1,1);
plot(f,abs(data_f));
title('Signal''s spectrum');
subtitle('Magnitude');
xlabel('f (Hz)');
ylabel('Signal''s spectrum magnitude (V)');
grid on;
subplot(2,1,2);
plot(f,angle(data_f));
subtitle('Phase');
xlabel('f (Hz)');
ylabel('Signal''s spectrum phase (rad)');
grid on;

Filter Design Parameters

Thus, we conclude that we need a Low Pass filter with a pass frequency (Fp) of 1000 Hz and a reject frequency (Fr) of 1500 Hz. That gives us a transition bandwidth of 500 Hz and a cut-off frequency (fc) of 1250 Hz.

fc = 1250 / fs;
wc = 2*pi*fc;
TBW = 500 / fs;

FIR Filter Design using Window Method

We choose the window we want and calculate the coefficients vector k for the window and the filter. This vector k must be between -"half" and +"half" the length K of the filter/window. Remember that vector k needs a central element with value '0'; for this, we include some additional calculations.

K = floor(5.5 / TBW);
K = K + 1 - mod(K,2);
k_center = (K-1)/2+1;
k = (-(K-1)/2:(K-1)/2);
W = 0.42 + 0.5*cos(2*pi*k/(K-1)) + 0.08*cos(4*pi*k/(K-1));

Generate Ideal Filter

We generate the "ideal" filter based on the sinc() function.

hd_lp = 2*fc*(sin(2*pi*fc*k)./(2*pi*fc*k));
hd_lp(k_center) = 2*fc;

Apply Window and Create Real Filter

Then we apply the window function, and we get our real filter. And we calculate the frequency response of the filter using the FFT.

h_lp = hd_lp .* W;
H = fft(h_lp,nFFT);

Represent Filter Response

Now we represent the filter's response graphically.

figure;
stem(k,h_lp);
title('Filter''s impulse response');
xlabel('k');
ylabel('h[k]');
grid on;

figure;
subplot(2,1,1);
plot(f,20*log10(abs(H)));
title('Filter''s frequency response');
subtitle('Magnitude');
xlabel('f (Hz)');
ylabel('|H(f)| (dB)');
grid on;
subplot(2,1,2);
plot(f,angle(H));
subtitle('Phase');
xlabel('f (Hz)');
ylabel('\angleH(f) (rad)');
grid on;

Filter Signal and View Output

Finally, we filter the input signal by means of the convolution, and we represent the output in both time and frequency domains.

y = conv(h_lp,data);
N = length(y);
n = 0:N-1;
nTs = n*Ts;
figure;
plot(nTs,y);
title('Output signal, TIME domain');
xlabel('nT_{S} (s)');
ylabel('y[nT_{S}] (V)');
grid on;

Y = fft(y,nFFT) / min(N,nFFT);
figure;
subplot(2,1,1);
plot(f,abs(Y));
title('Output''s spectrum');
subtitle('Magnitude');
xlabel('f (Hz)');
ylabel('Output''s spectrum magnitude (V)');
grid on;
subplot(2,1,2);
plot(f,angle(Y));
subtitle('Phase');
xlabel('f (Hz)');
ylabel('Output''s spectrum phase (rad)');
grid on;