Filter Design menggunakan Finite Impulse Response (FIR) Band Pass di MATLAB

 komponen frekuensi= 5hz, 7Hz, 9Hz

Amplitudo= 1 2 3

Sample interval, dt= 8msec dan panjang signal 2 second

FIR Low Pass (5hz)

Cut off fc= 6hz

Panjang filter= 13, 25, 49, 73

Membuat Filter Band Pass(FBP) dengan cut off fc1= 6hz dan fc2 8hz
Panjang filter= 13, 25, 49, 73

Aplikasikan filter yang anda design tersebut untuk mengambil hanya komponen signal dengan frekwensi f2= 7 Hz 

clear

clc

 

%Generate Signal

Fs= 1000; %Sampling Frequency

Ts= 1/Fs; %Sampling period or time step

dt= 0:Ts:2-Ts; %Signal Duration

f1= 5;

f2= 7;

f3= 9;

 

%Y= A Sin (2pi f t + theta)

Y1= 1*sin(2*pi*f1*dt);

Y2= 2*sin(2*pi*f2*dt);

Y3= 3*sin(2*pi*f3*dt);

 

Y4= Y1+Y2+Y3;

 

figure(1)

subplot(4,1,1)

plot(dt,Y1)

title('Signal 5Hz')

xlabel('Time')

ylabel('Amplitude')

subplot(4,1,2)

plot(dt,Y2,'g')

title('Signal 7Hz')

xlabel('Time')

ylabel('Amplitude')

subplot(4,1,3)

plot(dt,Y3,'b')

title('Signal 9Hz')

xlabel('Time')

ylabel('Amplitude')

subplot(4,1,4)

plot(dt,Y4,'r')

title('Composite Signal')

xlabel('Time')

ylabel('Amplitude')

%Fourier Transform

fs=1./dt;

N=length(t); %Sample

f=fs*(0:(N/2))/N;

yfft=fft(Y4).*dt;

 

%PLOT

cut=floor(length(f));

figure(2);

subplot(2,1,1);

plot(t,Y4,'r');

title('Composite Signal');

xlabel('Time');

ylabel('Amplitude');

subplot(2,1,2);

plot(f,abs(yfft(1:length(f))));

title('FFT Composite SIgnal');

xlabel('Frequency');

ylabel('Magnitude');

axis('tight');

xlim([0 cut]);

 

%Band Pass Filter

fh1=6; %frekuensi cut-off

fh2=8;

 

h13(1)=(2*fh2)-(2*fh1);

h25(1)=(2*fh2)-(2*fh1);

h49(1)=(2*fh2)-(2*fh1);

h73(1)=(2*fh2)-(2*fh1);

 

wh1=2*pi*fh1;

wh2=2*pi*fh2;

k1= 13; %Filter length

k2= 25;

k3= 49;

k4= 73;

 

h13(2:k1+1)=((2*fh2).*sin(wh2.*(1:k1).*dt)./(wh2.*(1:k1).*dt))-((2*fh1).*sin(wh1.*(1:k1).*dt)./(wh1.*(1:k1).*dt));

h25(2:k2+1)=((2*fh2).*sin(wh2.*(1:k2).*dt)./(wh2.*(1:k2).*dt))-((2*fh1).*sin(wh1.*(1:k2).*dt)./(wh1.*(1:k2).*dt));

h49(2:k3+1)=((2*fh2).*sin(wh2.*(1:k3).*dt)./(wh2.*(1:k3).*dt))-((2*fh1).*sin(wh1.*(1:k3).*dt)./(wh1.*(1:k3).*dt));

h73(2:k4+1)=((2*fh2).*sin(wh2.*(1:k4).*dt)./(wh2.*(1:k4).*dt))-((2*fh1).*sin(wh1.*(1:k4).*dt)./(wh1.*(1:k4).*dt));

 

ht13=dt*[h13(k1:-1:2) h13];

ht25=dt*[h25(k2:-1:2) h25];

ht49=dt*[h49(k3:-1:2) h49];

ht73=dt*[h73(k4:-1:2) h73];

 

dfr113=1/(length(ht13)*dt);

dfr125=1/(length(ht25)*dt);

dfr149=1/(length(ht49)*dt);

dfr173=1/(length(ht73)*dt);

 

fnyq=1/(2*dt);

fr113=(0:dfr113:fnyq);

fr125=(0:dfr125:fnyq);

fr149=(0:dfr149:fnyq);

fr173=(0:dfr173:fnyq);

 

cfr113=(1:length(fr113));

cfr125=(1:length(fr125));

cfr149=(1:length(fr149));

cfr173=(1:length(fr173));

 

fhlp13=fft(ht13);

fhlp25=fft(ht25);

fhlp49=fft(ht49);

fhlp73=fft(ht73);

 

fhlpc13=fhlp13(cfr113);

fhlpc25=fhlp25(cfr125);

fhlpc49=fhlp49(cfr149);

fhlpc73=fhlp73(cfr173);

 

%PLOT

figure(3);

subplot(2,2,1);

plot(fr113,abs(fhlpc13));

xlim([0 cut]);

title('BPF freq. domain k=13');

xlabel('Frequency');

ylabel('Magnitude');

subplot(2,2,2);

plot(fr125,abs(fhlpc25));

xlim([0 cut]);

title('BPF freq. domain k=25');

xlabel('Frequency');

ylabel('Magnitude');

subplot(2,2,3);

plot(fr149,abs(fhlpc49));

xlim([0 cut]);

title('BPF freq. domain k=49');

xlabel('Frequency');

ylabel('Magnitude');

subplot(2,2,4);

plot(fr173,abs(fhlpc73));

xlim([0 cut]);

title('BPF freq. domain k=73');

xlabel('Frequency');

ylabel('Magnitude');

 

kmin1=-1*k1;

kmin2=-1*k2;

kmin3=-1*k3;

kmin4=-1*k4;

tk13=kmin1+1:1:k1;

tk25=kmin2+1:1:k2;

tk49=kmin3+1:1:k3;

tk73=kmin4+1:1:k4;

 

%PLOT

figure(4);

subplot(2,2,1);

plot(tk13,ht13);

title('BPF time domain k=13');

xlabel('Time');

ylabel('Amplitude');

subplot(2,2,2);

plot(tk25,ht25);

title('BPF time domain k=25');

xlabel('Time');

ylabel('Amplitude');

subplot(2,2,3);

plot(tk49,ht49);

title('BPF time domain k=49');

xlabel('Time');

ylabel('Amplitude');

subplot(2,2,4);

plot(tk73,ht73);

title('BPF time domain k=73');

xlabel('Time');

ylabel('Amplitude');

% %apply to time series

ylpf13=conv(Y4,ht13);

ylpf25=conv(Y4,ht25);

ylpf49=conv(Y4,ht49);

ylpf73=conv(Y4,ht73);

 

yylpf13=ylpf13(k1:length(ylpf13)-k1);

yylpf25=ylpf25(k2:length(ylpf25)-k2);

yylpf49=ylpf49(k3:length(ylpf49)-k3);

yylpf73=ylpf73(k4:length(ylpf73)-k4);

 

figure(5)

subplot(4,1,1)

plot(t,Y4,'k');hold on

plot(t,yylpf13,'r');

title('Perbandingan sinyal asli dgn filter k=13 BPF')

xlabel('time')

ylabel('Amplitude');

hold off

subplot(4,1,2)

plot(t,Y4,'k');hold on

plot(t,yylpf25,'r');

title('Perbandingan sinyal asli dgn filter k=25 BPF');

xlabel('time');

ylabel('Amplitude');

hold off

subplot(4,1,3)

plot(t,Y4,'k');hold on

plot(t,yylpf49,'r');hold on

title('Perbandingan sinyal asli dgn filter k=49 BPF');

xlabel('time');

ylabel('Amplitude');

hold off

subplot(4,1,4)

plot(t,Y4,'k');hold on

plot(t,yylpf73,'r');

title('Perbandingan sinyal asli dgn filter k=73 BPF');

xlabel('time');

ylabel('Amplitude');

hold off