[selina]

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.1;
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=254;
for i=1:n
yy1(i)= b(4)*y(i+3)+ b(3)*y(i+2)+ b(2)*y(i+1) + b(1)*y(i)+a(4)*yy1(i+2)+a(3)*yy1(i+1) +a(2)*yy1(i);
end
%滤波
y=fft(y,N); %将信号做FFT变换
pyy=y.*conj(y); %做功率谱分析
f=(0N/2-1));
figure(1);
plot(f,pyy(1:N/2))

y=fft(yy1,N); % 将滤波后数据做功率谱分析
pyy=y.*conj(y);
f=(0N/2-1));
figure(2);
plot(f,pyy(1:N/2))


a =

1.0000 1.1795 0.9180 0.2742

??? Attempted to access yy1(3); index out of bounds because numel(yy1)=2.

Error in ==> iir at 25
yy1(i)= b(4)*y(i+3)+ b(3)*y(i+2)+ b(2)*y(i+1) + b(1)*y(i)+a(4)*yy1(i+2)+a(3)*yy1(i+1) +a(2)*yy1(i);