功率谱复现不收敛
问题描述：（以上分别为迭代10次 迭代12次级迭代14次的结果）功率谱复现理论实验阶段，采用的方法是：1．采用二阶振荡系统模拟目标系统；2．随机定义时域信号（加入一定的噪声）求得目标功率谱，即参考谱；3．频响函数辨识：利用参考谱作为系统驱动谱，得到幅值信息后，采用频域随机化得到频谱，通过ifft得到时域信号，进行时域随机化后，真随机信号输入到系统中，得到输出谱，进行均衡后，再次得到驱动谱，不停的进行迭代，比较误差；4．考核方法：参考谱正负3db以内。发现：1.功率谱分析采用的是matlab自带的PSD函数，发现fft后的频谱幅值与PSD后的幅值关系为：Aw_drv=9.74*sqrt(Gx);2.单纯的使用传统的均衡迭代的方法，得到的复现功率谱并不是很好，于是在此基础上改进，采用每次均衡后得到相应的频响函数参与误差迭代；3.不管是基于那种方法，共同的问题是：迭代过程并不收敛，有时收敛的很好，有时很发散，不知道原因何在；4.加入时域随机化得到真随机信号时，功率谱很大很大，也许是“随机抽头、加窗、重叠”做法有问题，采用方法是随机抽头，半正弦窗后延时1/2个周期重叠。 function [ ] = Testpsd()%TESTPSD Summary of this function goes here%用到的函数有：discretePSD()、stopby()、tmRand()%%%%得到参考功率谱fs=32;nfft=64;ts=2;t=0:1/fs:ts;GaussNos=wgn(1,length(t),5);x_spd=2*sin(2*pi*2*t+pi/2)+4*sin(2*pi*5*t+pi/4)+6*sin(2*pi*10*t+pi/3)+GaussNos;figure;plot(t,x_spd);Rxx=psd(x_spd,nfft,fs,hanning(nfft),0,0.99);Rxx=Rxx';index=0:nfft/2;f=index*fs/nfft;%%%%产生系统辨识的驱动信号：由参考谱产生(驱动信号未知)Aw_drv=9.74*sqrt(Rxx);fai=unifrnd(-pi,pi,1,nfft/2-1);Aw_drv(2:end-1)=Aw_drv(2:end-1).*exp(fai*i);R_drv=;x_drv=real(ifft(R_drv,nfft)); %%伪随机信号%%%%%%%%%频响函数估计w=8*2*pi;ksai=0.2;num=;den=;G=tf(num,den);%%%%%%%%%输出信号y_drv=lsim(G,x_drv,t(1:end-1));y_drv=y_drv';%%%%参考谱values=spcrv([;[ Rxx(1) Rxx Rxx(end)]],3,1000);figure;plot(values(1,:),values(2,:),'r') ;%%%%%%%%生成初始的驱动信号Gx0=Rxx;x_drv0=discretePSD(Gx0);y_drv0=lsim(G,x_drv0,t(1:end-1));Gy0=psd(y_drv0,nfft,fs,hanning(nfft),0,0.99);Gy0=Gy0';Hff=Gy0./Gx0;values=spcrv([;],3,1000);hold onplot(values(1,:),values(2,:)) ;Ef0=Rxx-Gy0;for m=1:20 %%%%%画图 Rrr=Rxx;Rr1=Rrr*0.70794;Rr2=Rrr*1.4125; values=spcrv([;[ Rrr(1) Rrr Rrr(end)]],3,1000); figure;plot(values(1,:),values(2,:),'r') ; values=spcrv([;[ Rr1(1) Rr1 Rr1(end)]],3,1000); hold on;plot(values(1,:),values(2,:),'r--') ; values=spcrv([;[ Rr2(1) Rr2 Rr2(end)]],3,1000);hold on;plot(values(1,:),values(2,:),'r--') ;%%%%%%%% Gx1=Gx0+Ef0./Hff; x_drv0=discretePSD(Gx1);%频域随机化% x_drv0=tmRand( x_drv0);%时域随机化 y_drv0=lsim(G,x_drv0,t(1:end-1)); Gy0=psd(y_drv0,nfft,fs,hanning(nfft),0,0.99); Gy0=Gy0'; %%%%%%%%%%%%%values=spcrv([;],3,1000); hold on plot(values(1,:),values(2,:)) ; Hff=Gy0./Gx1; Ef0=Rxx-Gy0; Gx0=Gx1;endend%%%%%%%%%%%%%%%%%%%%频域随机化方法如下：functionx_drv = discretePSD( Gx)%DISCRETEPSD Summary of this function goes here% 频域随机化nfft=64;Aw_drv=9.74*sqrt(Gx);fai=unifrnd(-pi,pi,1,nfft/2-1);Aw_drv(2:end-1)=Aw_drv(2:end-1).*exp(fai*i);R_drv=;x_drv=real(ifft(R_drv,nfft)); %%%%接下来进行时域随机化end {:{39}:} 表示很强大，我不懂帮顶一下吧
页:
[1]