功率谱复现不收敛

ohzking

迭代10次

迭代12次

迭代14次

问题描述：（以上分别为迭代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=[Aw_drv,conj(Aw_drv(end-1:-1:2))];

x_drv=real(ifft(R_drv,nfft));   %%伪随机信号

%%%%%%%%%频响函数估计

w=8*2*pi;ksai=0.2;

num=[w^2];den=[1 2*ksai*w w^2];

G=tf(num,den);

%%%%%%%%%输出信号

y_drv=lsim(G,x_drv,t(1:end-1));

y_drv=y_drv';

%%%%参考谱

values=spcrv([[f(1) f f(end)];[ 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([[f(1) f f(end)];[Gy0(1) Gy0 Gy0(end)]],3,1000);

hold on

plot(values(1,:),values(2,:)) ;

Ef0=Rxx-Gy0;

for m=1:20

    %%%%%画图

    Rrr=Rxx;Rr1=Rrr*0.70794;Rr2=Rrr*1.4125;

    values=spcrv([[f(1) f f(end)];[ Rrr(1) Rrr Rrr(end)]],3,1000);

    figure;plot(values(1,:),values(2,:),'r') ;

    values=spcrv([[f(1) f f(end)];[ Rr1(1) Rr1 Rr1(end)]],3,1000);

    hold on;plot(values(1,:),values(2,:),'r--') ;

    values=spcrv([[f(1) f f(end)];[ 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([[f(1) f f(end)];[Gy0(1) Gy0 Gy0(end)]],3,1000);

    hold on

    plot(values(1,:),values(2,:)) ;

    Hff=Gy0./Gx1;

    Ef0=Rxx-Gy0;

    Gx0=Gx1;

end

end

%%%%%%%%%%%%%%%%%%%%

频域随机化方法如下：

function  x_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=[Aw_drv,conj(Aw_drv(end-1:-1:2))];

x_drv=real(ifft(R_drv,nfft));

%%%%接下来进行时域随机化

end

璀璨星辰

lbtv

表示很强大，我不懂帮顶一下吧