iir_matlab.m

DSP系统入门与实践源代码

clear all;                    	     %清寄存器值
clf;                              			% 清屏
N=256;                            			 %数据点数
fs=600;                            			% 采样频率
dt=1/fs;
                       
for k=1:N;                    
   f1=300;                         			 %信号频率
   f2=100;                         			 %信号频率
   y(k)=sin(2*pi*f1*k*dt)+sin(2*pi*f2*k*dt)+sin(2*pi*(f1+20)*k*dt);   %产生信号
end

lp=200;                             			 %截止频率
wn1=2*lp/fs;                        			 %函数的参数
[z1,p1,k1] = CHEBY1(3,0.5,wn1,'high');      			 %滤波器的极零点表示
[B,A] = CHEBY1(3,0.5,wn1,'high');                      	% 滤波器的传递函数表示

yy1(1)=0;
yy1(2)=0;
yy1(3)=0;
b(1)=0.0580; b(2)=-0.1741; b(3)=0.1741;b(4)=-0.058;
a(1)=1; a(2)=1.1795; a(3)=0.918;a(4)=0.2742;

n=253;
for   i=1:n
     yy1(i+3)= b(4)*y(i)+ b(3)*y(i+1)+ b(2)*y(i+2)+ b(1)*y(i+3)-a(4)*yy1(i)-a(3)*yy1(i+1)-a(2)*yy1(i+2);           
end
%yy1=filter(B,A,y);                   	%滤波
y=fft(y,N);                         	%将信号做FFT变换
pyy=y.*conj(y);                     	%做功率谱分析
f=(0:(N/2-1));                                 %坐标转换
figure(1);
plot(2.33*f,pyy(1:N/2))
                
y=fft(yy1,N);                      	%将滤波后数据做功率谱分析
pyy=y.*conj(y);                                %做功率谱分析
f=(0:(N/2-1));                                   %坐标转换
figure(2);
plot(2.33*f,pyy(1:N/2));                          %画功率谱