声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1482|回复: 2

功率谱复现不收敛

[复制链接]
发表于 2013-9-2 16:28 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x

迭代10次

迭代10次

迭代12次

迭代12次

迭代14次

迭代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
回复
分享到:

使用道具 举报

发表于 2014-6-23 15:26 | 显示全部楼层
发表于 2014-6-24 10:28 | 显示全部楼层
表示很强大,我不懂帮顶一下吧
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-24 20:30 , Processed in 0.087558 second(s), 22 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表