现代信号处理的大作业:LD算法以及WV变换 下载本文

noise = randn(1,len);% 产生均值为0,方差(功率)为1,数据长度为3000的高斯白噪声

a0=[1 0.72 0.88];% 已知模型的输入系数 x=filter(1,a0,noise);%

x(n)=a0(1)*noise(n-2)+a0(2)*noise(n-1)+a0(3)*noise(n);

%% 利用L-D算法预测模型参数 m = 50;

[a,P] = Levinson_Durbin(x,m); aa = [ones(m,1),a];% 预测出来的系数 %预测误差功率图 figure

plot(1:m,P(1:m),'rh') ylim([0,6])

title('预测误差功率随迭代次数的变化'); xlabel('迭代次数') ylabel('预测误差功率') % 利用aryule验证 y = aryule(x,m); figure subplot(1,2,1) plot(1:51,aa(m,:),'r')

title('编写程序运行得出的系数') subplot(1,2,2) plot(1:51,y,'k')

title('aryule运行得出的系数')

6

4.2、LD算法子函数: %% Levinson-Durbin algorithm function [a,P]=Levinson_Durbin(x,m) % x为输入信号 % m为阶数 a = zeros(m,m); N = length(x); P = zeros(1,m + 1); for i = 1:N

P(1) = P(1) + x(i)^2; end

P(1) = P(1)/N;% 计算预测误差功率的初始值P_0

f = zeros(m + 1,N);% m阶前向预测误差 g = zeros(m + 1,N);% m阶后向预测误差 f(1,:) = x; g(1,:) = x; K = zeros(1,m); % 开始迭代: for i = 1:m numerator = 0; denominator = 0; for j = i + 1:N

numerator = numerator - f(i,j)*g(i,j - 1);

denominator = denominator + (f(i,j)^2 + g(i,j - 1)^2)/2; end

K(i) = numerator/denominator;% 反射系数 for j = 1:i - 1

7

a(i,j) = a(i - 1,j) + K(i)*a(i - 1,i - j); end a(i,i) = K(i); if i >= 2

P(i) = (1 - K(i)^2) * P(i - 1);% 更新预测误差功率 end for j = i + 1:N

f(i + 1,j) = f(i,j) + K(i)*g(i,j-1); g(i + 1,j) = K(i)*f(i,j) + g(i,j-1); end end

4.3、仿真结果:

显然,我们想要预测估计的是a0= [1 0.72 0.88],预测结果为(预测了1-50阶,后面阶数过多不易展示,因此只列出前1-6阶预测结果):

表1:LD算法预测系数表(1-6阶)

1 1 1 1 1 1 0.382709 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.719284 0.879451 0.688918 0.854616 -0.03453 0.688855 0.856167 -0.03328 0.001816 0.688871 0.855869 -0.02561 0.007983 0.008953 0.689077 0.856053 -0.0262 0.027638 0.024773 0.022965 以及预测误差功率随阶数变化图:

8

图1 、LD算法预测误差功率随阶数的变化

显然从2阶开始预测误差功率就显著减小几乎不变。用自带函数aryule计算系数,并对比50阶的各系数(横坐标为第几个系数,纵坐标为系数的值):

图2 、预测系数的对比

可以验证LD算法程序是有效的。

9

二.WV变换

1、问题描述

建立一个非平稳信号,由三个线性调频信号叠加而成:

?1???z(t)?()4[exp(?(t?t1)2?j?1t)?exp(?(t?t2)2?j?2t)?exp(?(t?t3)2?j?3t)]?222其中t1?t2?t3,?1??2??3。

求z(t)的WV分布,指出并分析其信号项和交叉项。

2、WV分布的分析

信号z(t)的Winger-Ville分布定义为:

W(t,w)??11z(t??)z*(t??)e?j?wd? (1) ??22?由二次叠加原理,Wigner-Ville分布的信号项和交叉项分别为:

??Wauto(t,w)?Wz1(t,w)?Wz2(t,w)???Wcross(t,w)?Wz1,z2(t,w)?Wz2,z1(t,w) (2)

其中,

1/4?????2?z1???exp[?(t?t1)?jw1t]2???? (3) ?1/4?????2z?exp[?(t?t)?jw2t]2??1?2????

由WV分布的定义得:

*Wz1(t,w)??z1(t?)z1(t?)e?j?wd???22????()1/2?exp[??2?j(w?w1)???(t?t1)2]d????4 (4)

???在积分公式

10