matlab使用invfreqs出错,数字信号实验报告(matlab实验数据与结果)

可见，u越大，h(n)和s(n)的速度越快。 4.2 自适应AR模型

输入估计AR模型的阶数M =2 输入实际AR模型的阶数P= 2 输入步长u =0.01 输入 W的方差1 输入样本个数 L= 100 输入a(1) =-1.3 输入a(2) =0.8

可见a1(n)、a2(n)随着n的增加趋近于a1、a2。

附录

实验源程序

⑴%用脉冲响应不变法设计IIR滤波器 clear all;close all;clc;

T=0.25;Fs=4;wc=2*pi;Omegac=wc/T;N=4;t=[0:0.01:3];n=[0:1:3/T]; [b,a]=butter(N,Omegac,'s'); [c,d]=impinvar(b,a,T); [h,w]=freqs(b,a); [ha,x,t]=impulse(b,a); mag = abs(h);mag = 20*log10(mag); figure(1);impz(ha);

title('h(n)');ylabel('幅度');grid on; figure(2)

plot(t,ha);title('单位脉冲响应ha(t)');grid on; figure(3)

subplot(1,2,1);plot(w,mag);title('幅度响应Ha(j/Omega)');grid on; [h,w]=freqz(b,a);mag = abs(h);mag = 20*log10(mag);

subplot(1,2,2);plot(w,mag);title('幅度响应H(e^jw)');grid on; %用Butterworth设计

Fs=1;Omegap=0.4;Omegar=0.6;Ap=1;Ar=25;T=1/Fs; [n,Wn] = buttord(Omegap,Omegar,Ap,Ar)

[cs,ds]= butter(n,Wn); [b,a]=impinvar(cs,ds,T); [h,w]=freqz(b,a); mag = abs(h);db = 20*log10(mag);

subplot(1,2,1);plot(w/pi,mag);title('幅度响应H(e^jw)');grid on;

subplot(1,2,2);plot(w/pi,db);title('巴特沃思型幅度特性(db)');grid on; %用chebyshev设计

Fs=1;Omegap=0.4*pi;Omegar=0.6*pi;Ap=1;Ar=25;T=1/Fs; [n,Wn] = cheb1ord(Omegap/pi,Omegar/pi,Ap,Ar); [cs,ds] = cheby1(n,Ap,Wn);

[b,a]=impinvar(cs,ds,T); [h,w]=freqz(b,a); mag = abs(h);db = 20*log10(mag);

subplot(1,2,1);plot(w/pi,mag);title('幅度响应H(e^jw)');grid on;

subplot(1,2,2);plot(w/pi,db);title('切比雪夫1型幅度特性(db)');grid on;

⑵%维纳滤波器 clear

L=input('输入样本点个数 ');N=input('输入FIR滤波器阶数 ');K=50;n=1:L; while(1)

w=sqrt(1-0.95^2)*(randn(L,1)); %获得高斯白噪声w

s(1)=w(1);v=randn(L,1); %获得方差为1的高斯白噪声v for i=2:L

s(i)=0.95*s(i-1)+w(i); %获得有用信号s(n) .公式s(n) = a*s(n-1)+w(n) end

for i=1:L

x(i) = s(i)+v(i); %随机信号x(n)

end; SI(1)=0; for i=2:L

SI(i)=0.724*SI(i-1)+0.238*x(i); end

%%%%%估计相关函数r_xx和r_xs%%%%%%%%%%%%%%% for i=1:K+1

rxx(i)=sum(x(i:L).*x(1:L-i+1))/(L-i+1); rxs(i)=sum(x(i:L).*s(1:L-i+1))/(L-i+1); end

r_xx_g=[rxx(K+1:-1:2),rxx(1:K+1)];r_xs_g=[rxs(K+1:-1:2),rxs(1:K+1)]; %%%%检验x的r_xx和r_xs是否与理论值相符%%%%%% r_xx_t = 0.95.^abs([-K:K]);r_xx_t(K+1) = r_xx_t(K+1)+1; r_xs_t = 0.95.^abs([-K:K]);

rou_xx=(sum((r_xx_g-r_xx_t).^2))/sum(r_xx_t.^2); rou_xs=(sum(r_xs_g-r_xs_t).^2)/sum(r_xs_t.^2); if rou_xx<0.03 & rou_xs<0.01 break; end end

n = 0:N-1;hi= 0.238*(0.724).^n; %求得理想的h(n) for i=1:N

xx(i)=sum(x(i:L).*x(1:L-i+1))/(L-i+1); end

for i=1:N

Rxx(i,1:N)=[xx(i:-1:1),xx(2:N+1-i)]; end

for i=1:N

xs(i)=sum(x(i:L).*s(1:L-i+1))/(L-i+1); end

invRxx=inv(Rxx);hr=invRxx*xs';SR=conv(hr,x); subplot(1,2,1)

plot(L-100:L,s(L-100:L),L-100:L,x(L-100:L),':'); title('s(n) —— x(n) ...'); figure(1)

subplot(1,2,2)

plot(L-100:L,s(L-100:L),L-100:L,SI(L-100:L),':'); title('s(n) —— SI(n) ...'); figure(2)

subplot(1,2,1)

plot(L-100:L,s(L-100:L),L-100:L,SR(L-100:L),':'); title('s(n) —— SR(n) ...'); subplot(1,2,2)

plot(1:N,hi(1:N),1:N,hr(1:N ),':');