ur filter

clear all

NN = 200;
dt = 0.0005;
TT = (0:dt:dt*(NN-1));

f1 = 50;
f2 = 100;
f3 = 150;
f4 = 200;
x = zeros(1,NN);
for n=1:NN
x(n) = 0.5*sin(2*pi*f1*dt*n) ...
     + 0.5*sin(2*pi*f2*dt*n) ...
     + 0.5*sin(2*pi*f3*dt*n) ...
     + 0.5*sin(2*pi*f4*dt*n);
end
figure(1);
subplot(3,2,1)
plot(TT,x,'k')
xlabel('T (sec)')
ylabel('x[k]')
axis( [ 0 dt*(NN-1)  -1.5 1.5 ])
set(gca,'fontsize',12)
grid on
title('FIR Filter')

% ------ Take the FFT

del_f = 1/(dt*NN);
FF = (0:del_f:del_f*(NN-1));

XX = (1/NN)*fft(x);
subplot(3,2,2)
plot(FF,abs(XX),'k')
xlabel('Freq (Hz)')
ylabel('X(z)')
axis( [ 0 250  0 .5 ])
set(gca,'XTick',[ 0 f1 f2 f3 f4 ])
set(gca,'fontsize',12)

% Frequency domain filter

%  Hanning
HH = zeros(1,NN);

HH = zeros(1,NN);
for n=1:6
HH(n) = 0.;
HH(NN-n+1) = 0.;
end

MM = 30;
for n=3:34
HH(n) = 0.5*(1 - cos(2*pi*(n-5)/MM));
HH(NN-n+1) = HH(n);
end

for n=49:NN/2
HH(n) = 0.;
HH(NN-n+1) = HH(n); % Notice that this supplies the negative frequencies.
end

% --- Mutliply XX by HH
YY = zeros(1,NN);
for n=1:NN
YY(n) = HH(n)*XX(n);
end

subplot(3,2,3)
plot(FF,HH,'k--')
hold on
plot(FF,abs(YY),'k');
axis( [ 00 250  0 1.1 ])
hold off
xlabel('Freq (Hz)')
ylabel('H(z),X(z)')
set(gca,'fontsize',12)
%legend('H','Y')

% --- Inv FFT and plot

y = NN*real(ifft(YY));
subplot(3,2,4)
plot(TT,y,'k')
ylabel('y[k]')
xlabel('T (sec)')
set(gca,'fontsize',12)

% Get the FIR filter from inverse HH
h = ifft(HH);

subplot(3,2,5)
plot(real(h),'ok');
axis( [ 1 8 -.2 .2  ])
grid on
ylabel('h[k]')
xlabel('k')
set(gca,'fontsize',12)

subplot(3,2,6)
plot(real(h),'ok');
axis( [  193 200 -.2 .2 ])
grid on
ylabel('h[k]')
xlabel('k')
set(gca,'fontsize',12)

saveas(gca,'fft.png')

figure(2);

% ---------------------------------------

% This part tests the FIR filter

% First, use h to develop a small FIR filter
hh = zeros(1,17);
hh(8) = real(h(1));
for n=2:7
hh(7+n) = real(h(n));
end

for n=1:8
hh(n) = hh(16-n);
end

subplot(3,2,1)
plot(hh,'ok')
axis( [ 1 15 -.3 .3 ])
grid on
xlabel('k')
ylabel('h[k]')

zz = conv(x,hh);       % Convolve the input with hh

for n=1:NN
z(n) = zz(n);
end

subplot(3,2,3)
plot(TT,z,'k')
axis( [ 0 .2 -1 1 ])
grid on
xlabel('sec')
ylabel('y[k]')

% Plot the FFT of the output
subplot(3,2,5)
plot(FF,(1/NN)*abs(fft(z)),'k');
axis( [ 0 200 0 .3 ])
set(gca,'XTick',[ 0 f1 f2 f3 ])
xlabel('Freq (Hz)')
ylabel('Y(z)')