现代信号处理的大作业LD算法以及WV变换.docx

1、现代信号处理的大作业LD算法以及WV变换现代信号处理大作业姓名： 学号：专业：电子科学与技术一. L-D算法的仿真实现1、问题描述用Matlab实现Levinsion-Durbin算法。2、算法分析2.1 、LD算法原理由于语音样点之间存在相关性，所以可以用过去的样点值来预测现在或未来的样点值。如下图所示图1 线性预测图示由上图可得，从而可以通过使实际语音x（n）和线性预测结果之间的误差e（n）在某个准则下达到最小值来决定唯一的一组预测系数。而这组系数就能反映语音信号的特性，可以作为语音信号特征参数来用于语音编码、语音合成和语音识别等应用中去。由估计值和实际信号值的误差,可有根据e(n)最小均

2、方误差准则，来决定唯一的一组预测系数，即:,由此可得到Y-W方程：，取遍k值之后有以下：由相关函数的偶函数性质有：在已知自相关函数的前提下，根据e(n)均方误差最小的原则来求解a，本实验中采用Levinson-Durbin算法。2.2、 LD算法的实现Levinson-Durbin算法首先由一阶AR模型开始，按照前面的Y-W方程可有，一阶AR模型（p=1）的Y-W方程是然后增加一阶，即令p=2，可得到：由上式可解出：然后令以此类推，可以得到一般的递推公式：式中的称为反射系数，而和是预测误差的均方误差值，因此必须大于等于0，这样应满足，进而得到，即预测误差随递推次数增加而减少。3、程序实现思路3

3、.1、观测值的选取。将均值为0.方差为1，数据长度3000的白噪声信号通过线性系统：H(z)=1/(1+a1*z(-1)+a2*z(-2)，设置模型的输入系数为a0= 1 0.72 0.88；由输出信号作为观测信号x(n)。3.2、通过LD算法函数估计从一阶到50阶的系数。3.3、用Matlab自带的aryule函数预测50阶系数，与自己编写的程序结果对比。4、程序如下：4.1、主函数：clc;clear;close all;% 给出一个已知的模型，让已知信号经过该模型之后利用函数估计模型系数% 已知模型设为x(n)=a1*noise(n-2)+a2*noise(n-1)+a3*noise(n

4、)len = 3000;% 数据长度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;% 预测出来的系数%预测误差功率图figureplot(1:m,P(1:m),rh)ylim(0,6)

5、title(预测误差功率随迭代次数的变化);xlabel(迭代次数)ylabel(预测误差功率)% 利用aryule验证y = aryule(x,m);figuresubplot(1,2,1)plot(1:51,aa(m,:),r)title(编写程序运行得出的系数)subplot(1,2,2)plot(1:51,y,k)title(aryule运行得出的系数)4.2、LD算法子函数：% Levinson-Durbin algorithmfunction a,P=Levinson_Durbin(x,m)% x为输入信号% m为阶数a = zeros(m,m);N = length(x);P =

6、 zeros(1,m + 1);for i = 1:N P(1) = P(1) + x(i)2;endP(1) = P(1)/N;% 计算预测误差功率的初始值P_0f = 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 =

7、 denominator + (f(i,j)2 + g(i,j - 1)2)/2; end K(i) = numerator/denominator;% 反射系数 for j = 1:i - 1 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-

8、1); endend4.3、仿真结果：显然，我们想要预测估计的是a0= 1 0.72 0.88，预测结果为（预测了1-50阶，后面阶数过多不易展示，因此只列出前1-6阶预测结果）：表1：LD算法预测系数表（1-6阶）10.3827090000010.7192840.879451000010.6889180.854616-0.0345300010.6888550.856167-0.033280.0018160010.6888710.855869-0.025610.0079830.008953010.6890770.856053-0.02620.0276380.0247730.022965以及预测

9、误差功率随阶数变化图：图1 、LD算法预测误差功率随阶数的变化显然从2阶开始预测误差功率就显著减小几乎不变。用自带函数aryule计算系数，并对比50阶的各系数（横坐标为第几个系数，纵坐标为系数的值）：图2 、预测系数的对比可以验证LD算法程序是有效的。二WV变换1、问题描述建立一个非平稳信号，由三个线性调频信号叠加而成：其中。求的WV分布，指出并分析其信号项和交叉项。2、WV分布的分析信号z(t)的Winger-Ville分布定义为： (1)由二次叠加原理，Wigner-Ville分布的信号项和交叉项分别为： (2)其中， (3)由WV分布的定义得： (4)在积分公式 (5)中，令和，得：

10、(6)同理，有： (7)故WV分布的信号项为： (8)又由交叉项的定义知： (9)利用积分公式，代入相应的参数，得： (10)故WV分布的交叉项为： (11)式中 (12)则信号的WV分布为： (13)3、程序设计close all;format long;syms x t;aifa=10;%设置初始参数aifa=input(Please input aifa:); %t1=input(Please input t1:);t2=input(Please input t2:);t3=input(Please input t3(t1t2t3):);w1=input(Please input w1:

11、);w2=input(Please input w2:); w3=input(Please input w3(w1w2w3):); %求取z(t+x/2)z1=exp(-aifa/2*(t+x/2-t1)2+1j*w1*(t+x/2)+exp(-aifa/2*(t+x/2-t2)2+1j*w2*(t+x/2)+exp(-aifa/2*(t+x/2-t3)2+1j*w3*(t+x/2); %求取z(t-x/2)的共轭 z2=exp(-aifa/2*(t-x/2-t1)2-1j*w1*(t-x/2)+exp(-aifa/2*(t-x/2-t2)2-1j*w2*(t-x/2)+exp(-aifa/2

12、*(t-x/2-t3)2-1j*w3*(t-x/2); %为求WV分布信号项作准备z3=exp(-aifa/2*(t+x/2-t1)2+1j*w1*(t+x/2)*exp(-aifa/2*(t-x/2-t1)2-1j*w1*(t-x/2); z4=exp(-aifa/2*(t+x/2-t2)2+1j*w2*(t+x/2)*exp(-aifa/2*(t-x/2-t2)2-1j*w2*(t-x/2);z5=exp(-aifa/2*(t+x/2-t3)2+1j*w3*(t+x/2)*exp(-aifa/2*(t-x/2-t3)2-1j*w3*(t-x/2);z=z1*z2;Z_auto=z3+z4+

13、z5;Wz=simplify(fourier(z)*(aifa/pi)(1/4)Wz_auto=simplify(fourier(Z_auto)*(aifa/pi)(1/4)Wz_cross=Wz-Wz_autoWz=inline(Wz)Wz_auto=inline(Wz_auto)t=0:0.01:3;f=-3.0:0.01:3;T,F=meshgrid(t,f);W=double(Wz(T,2*pi*F);W_auto=double(Wz_auto(T,2*pi*F); figure(name,WV分布);mesh(T,F,W); figure(name,WV分布信号项);mesh(T,F,W_auto); W_cross=W-W_auto;figure(name,WV分布交叉项);mesh(T,F,W_cross);4、程序运行结果1）值时，WV分布为图1 、WV分布图3 、WV分布的信号项（左）与交叉项（右）2. 值时，WV分布为图4、WV分布图5、WV分布的信号项（左）与交叉项（右）