马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
问题描述:(以上分别为迭代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 |