example5_4.m

《高阶谱分析》例题5.4

% Example 5.4
clc
clear

% 系统函数
b1=[1,-0.9];b2=[1,-1.52/0.617,1/0.617];
b=conv(b1,b2);  % 构建零点
a1=[1,-0.86];a2=[1,-1,0.61];
a=conv(a1,a2);  % 构建极点
figure(1)
zplane(b,a);    % 根据零极点画出Z平面
xlabel('Re');ylabel('Im');

% 系统函数的频率响应 
[H,W]=freqz(b,a,128);
H=H.*(0.617*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=0.617*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+1)+B(1-j)*out(k+q+1-j+1);     % 右移q+1位
    end
    out(k+q+1)=out(k+q+1)/k;
end

% 估计系统的频率相应
clear j
H_i=zeros(1,128);
H_o=zeros(1,128);
for i=1:128
    omega=(i-1)*pi/128;
    temp=exp(j*omega);
    for k=0:p
        H_i(i)=H_i(i)+in(k+1)*temp^(-k);
    end
    for k=-q:0
        H_o(i)=H_o(i)+out(k+q+1)*temp^(-k);
    end
end
H_sys=H_i.*H_o;
mag=abs(H_sys);
ang=angle(H_sys);
% end

figure(2)
hold on;
plot(W,mag,'r');
xlabel('\omega');ylabel('|H(\omega)|')
figure(3)
hold on;
plot(W,ang,'r');
xlabel('\omega');ylabel('\Phi_{h}(\omega)')
figure(4)
stem(1:p,A,'.');
xlabel('k');ylabel('A^{(k)}');
figure(5)
stem(1:q,B,'.');
xlabel('k');ylabel('B^{(k)}');