Matlab FIR, IIR fiter 설계 및 simulation 예제

교육 4일차에 직접 짠 Matlab code 이다.
아래 code는 두개의 wave file 을 읽어서  벌 소리를 filtering 하는 filter 를 만들고 ( 계수를 뽑아낸다는) 그래서 filtering 하고 fft 를 통해서 frequency domain 을 관찰하고 sound 로 들어보는 
matlab code이다.
%

clear all;
close all;


fs = 10000;
fs1 = 10000;
fnm = fs/2;
M = 40 ;
wn = 1000/fnm;
%b =  fir1(M, fnm);
b = fir1(M,wn);
b'

junhyuk_board_filtering.m

clear all
close all 

[x1, fs1] = wavread('../사운드 필터링 실습 화일/bul_1.wav');
[x2, fs2] = wavread('../사운드 필터링 실습 화일/bulpokpung_1.wav');


M = 50

Ts1 = 1/fs1;
fnm1 = fs1/2;
N1 = length(x1);
t1 = 0:Ts1:(N1-1)*Ts1;

subplot(5,1,1), plot(t1,x1)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%bul_1.wav%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
hertz1=(t1/Ts1)*(fs1/N1);

Res = fft(x1); magRes = abs(Res);
figure(2)
subplot(2,1,1), plot(hertz1(1:N1/3), magRes(1:N1/3))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%bulpokpung_1.wav%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ts2 = 1/fs2;
fnm2 = fs2/2;
N2 = length(x2);
t2 = 0:Ts2:(N2-1)*Ts2;
hertz2=(t2/Ts2)*(fs2/N2);

Res = fft(x2); magRes = abs(Res);
subplot(2,1,2), plot(hertz2(1:N2/3), magRes(1:N2/3)) %x time freq plot %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%bulpokpung_1.wav%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


f = [0 100/fnm2 100/fnm2   400/fnm2 400/fnm2 2800/fnm2 2800/fnm2 1];
a = [0  0       0        0.5      1       1 1 1];
b= fir2(M,f,a);
b'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%high pass filtering%%%%%%%%%%%%%%%%%%%%%%%%%%
result = filter(b,1,x2);
Res = fft(result); magRes = abs(Res);
figure(4)
subplot(2,1,1), plot(hertz1(1:N2/2), magRes(1:N2/2)) %x time freq plot %
soundsc(result,fs2)
pause 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%Low pass filtering%%%%%%%%%%%%%%%%%%%%%%%%%%
f = [0 500/fnm2 500/fnm2 1000/fnm2 1000/fnm2 1];


a = [1  1       1        0      0       0];
b= fir2(M,f,a);

result = filter(b,1,x2);
Res = fft(result); magRes = abs(Res);
figure(4)
subplot(2,1,2), plot(hertz1(1:N2/2), magRes(1:N2/2)) %x time freq plot %
soundsc(result,fs2)
pause