,求冲激响应。.doc

高阶谱分析2007年试题
姓名:薛续磊学号:SA06006104 成绩:
。
例题的计算流程如下:
构建题目中的系统传递函数,求冲激响应。
产生零均值非高斯白噪声作为系统激励,并计算系统输出。
按照下式计算输出的三阶矩:
()
其中
()
构建如下方程
()
其中为的矩阵,其元素形如; ,大小为;是的矢量,其元素形如。
然后使用最小二乘法计算:
()
则由上式所求得的可得到倒谱参数和,则如下计算和
()
()
由和可得和,则
()
程序源代码如下:
clc
clear
% 系统函数
b1=[1,-];b2=[1,-,1/];
b=conv(b1,b2); % 构建零点
a1=[1,-];a2=[1,-1,];
a=conv(a1,a2); % 构建极点
figure(1), zplane(b,a); % 根据零极点画出Z平面
xlabel('Re');ylabel('Im');
% 系统函数的频率响应
[H,W]=freqz(b,a,128);
H=H.*(*exp(j*W).^2);
figure(2), plot(W,abs(H));
figure(3), plot(W,angle(H));
% 产生系统激励
% for num=1:100 % 100次平均
K=32; % 数据段的数目
M=2048;
for i=1:K*M
t=rand;
if t<=1/3
x(i)=2*4^(1/3)*rand;
else
x(i)=-4^(1/3)*rand;
end
end
y=filter(b,a,x);
y=y(1:K*M); % 系统输出
y=*y;
% 估计系统输出的三阶统计量
p=20; q=20;
w=max(p,q); z=floor(w/2);
m_3=zeros(4*w+1,2*z+p+q+1);
for i=1:K % 共K段记录
m_temp=zeros(4*w+1,2*z+p+q+1);
y0=y((i-1)*M+1:i*M); % 取M长的一段
for n=-2*w:2*w
for l=-z-q:z+p
s1=max([0,-n,-l]);
s2=min([M-1,M-1-n,M-1-l]);
for j=s1:s2
temp=y0(j+1)*y0(j+1+n)*y0(j+1+l);
m_temp(n+2*w+1,l+z+q+1)=m_temp(n+2*w+1,l+z+q+1)+temp;
end
m_temp(n+2*w+1,l+z+q+1)=m_temp(n+2*w+1,l+z+q+1)/M; % 该段记录的3阶矩
end
end
m_3=m_3+m_temp; % K段累加
end
m_3=m_3/K; % K段平均
% 计算倒谱参数
U=zeros((2*w+1)*(2*z+1),(p+q));
u=zeros((2*w+1)*(2*z+1),1);
count=0;
for n=-w:w
for l=-z:z
count=count+1;
u(count,1)=-n*m_3(n+2*w+1,l+z+q+1);
for i=1:p+q
if i<=p
U(count,i)=m_3(n-i+2*w+1,l+z+q+1)-m_3(n+i+2*w+1,l+i+z+q+1);
else
temp=i-p;
U(count,i)=m_3(n-temp+2*w+1,l-temp+z+q+1)-m_3(n+temp+2*w+1,l+z+q+1);
end
end
end
end
v=inv(U'*U)*U'*u;
for i=1:p+q
if i<=p
A(i)=v(i);
else
temp=i-p;
B(temp)=v(i);
end
end
% 计算i(k)和o(k)
in=zeros(1,p+1); out=zeros(1,q+1);
in(1)=1; out(q+1)=1;
for k=1:p
for j=2:k+1
in(k+1)=in(k+1)+A(j-1)*in(k+1-j+1); % 右移1位
end
in(k+1)=-1/k*in(k+1);
end
for k=-1:-1:-q
for j=k+1:0
out(k+q+1)=out(k+q