正弦扫频信号幅值及相位的提取

1正弦扫频信号幅值及相位的提取(1)正弦振动控制系统提供输入的扫频信号，对于对数扫频，，其中Sr为对数扫描率，若频响函数为则系统输出为。测量系统中可得到Calo信号及响应信号，通过对二者进行数据处理，可得到频域下的响应。不知道LMS的信号采集软件是如何提取频域响应的，个人认为软件计算速度有限，LMS应该是通过硬件实现的。下面我提供几种方法并进行比较。算例对于Calo信号，频响函数为，其中，信号采样率为1000次/秒，图1给出了时域下的响应信号。图1时域下的响应信号正弦扫频信号幅值及相位的提取(2)2方法1分段FFT在[f,f+df]区间内对Calo信号、响应信号进行FFT变换，二者在频率f处的谱值比即为频响函数在f处的值。此方法的缺陷是由于信号采样率为1000Hz，而[f,f+df]的区间很窄，在此区间下时域的点不会很多，因而FFT的频率分辨率不高。对于没有相位差的扫频信号，此方法能较好的提取幅值。图2给出了使用此方法提取的幅值与理论结果比较，由图中可以看出二者基本吻合。图2使用分段FFT提取的频域幅值对于有相位差的扫频信号，则要对结果进行光滑处理，Matlab的smooth函数提供了这一功能。图3给出了有相位差时分段FFT提取的幅值与相位同理论结果的比较，从图中可以看出在频域峰值处分段FFT比理论值大，在其余频段二者吻合较好。3图3使用分段FFT提取的频域幅值、相位下面给出了实现分段FFT提取扫频信号的频域幅值、相位的Matlab代码。----------------------------------------------------------------------%Decomposetheamplitudeandphasefromthesweepsignal%Localfftandsmoothareemployed.f1=5;%theinitialfreqs=4;%sweepratefr=50;%Resonantfreqaf=[];%amplitudepf=[];%phasek1=0.02;%dampingratiodf=0.01;%freqintervalforfa=40:df:60t1=60/s*log2(fa/f1);t2=60/s*log2((fa+df)/f1);ta=t1:0.001:t2;N=length(ta);4ft=f1*2.^(s/60*ta);A1=sin(2*pi*ft.*ta);lamb=ft/fr;B1=1./(1-lamb.^2+j*2*k1*lamb);%transferfunctionA2=abs(B1).*sin(2*pi*ft.*ta+angle(B1));ffreq=exp(-j*2*pi*(fa-400)*ta);%freqshiftfortimedomainspa=fft(ffreq.*A1);spb=fft(ffreq.*A2);spr=abs(spb./spa);spp=angle(spb./spa);k=ceil(N*0.001*400);af=[af,spr(k+1)];pf=[pf,spp(k+1)];endaf=smooth(af,7);pf=smooth(pf,7);%Keytrickfa=40:df:60;lamb=fa/fr;bf=abs(1./(1-lamb.^2+j*2*k1*lamb));subplot(2,1,1);plot(fa,af,'r-',fa,bf,'b-.');legend('NumericResult','TheoreticResult');title('AmplitudeofSweepSignal');xlabel('f');ylabel('A(f)');subplot(2,1,2);bpf=angle(1./(1-lamb.^2+j*2*k1*lamb));plot(fa,180/pi*pf,'r-',fa,180/pi*bpf,'b-.');legend('NumericResult','TheoreticResult');title('PhaseofSweepSignal');5xlabel('f');ylabel('\Psi(f)');-----------------------------------------------------------------------分段FFT提取方法计算速度一般，不会出现异常而中止，计算精度基本也能保证。正弦扫频信号幅值及相位的提取(3)方法2分段曲线拟合在[f,f+df]区间内，假定A，ψ不变，此区间内在时域内对其拟合。图4给出了有相位差时曲线拟合提取的幅值与相位同理论结果的比较，从图中可以看出计算结果与真实值吻合非常好。图4使用分段曲线拟合提取的频域幅值、相位分段曲线拟合提取的结果精度非常高，但是由于是拟合方法，因而可能会由于初始值给的不合理或拟合关系式不恰当而出现迭代次数超过规定值从而导致计算中止。下面给出了实现分段曲线拟合提取扫频信号的频域幅值、相位的Matlab代码。---------------------------------------------------------------6%Decomposetheamplitudeandphasefromthesweepsignal%Localcurvefitisapplied.f1=5;%theinitialfreqs=4;%sweepratefr=50;%Resonantfreqaf=[];%amplitudek1=0.02;%dampingratiodf=0.01;%freqincreasea=[];%amplitudeph=[];%phasex0=[1,0];%initialguessforfa=40:df:60;t1=60/s*log2(fa/f1);t2=60/s*log2((fa+df)/f1);ta=t1:0.001:t2;N=length(ta);ft=f1*2.^(s/60*ta);lamb=ft/fr;B1=1./(1-lamb.^2+j*2*k1*lamb);%transferfunctionA1=abs(B1).*sin(2*pi*ft.*ta+angle(B1));xd=2*pi*ft.*ta;opt=optimset('Display','off');x=lsqnonlin('fsin',x0,[0,-pi],[inf,pi],opt,xd,A1);x0=x;%Keytricka=[a,x(1)];ph=[ph,x(2)];endft=40:df:60;lamb=ft/fr;B1=1./(1-lamb.^2+j*2*k1*lamb);%transferfunction7subplot(2,1,1);plot(ft,a,'r-',ft,abs(B1),'b-.');title('MagnitudeofSweepSignal');legend('NumericResult','TheoreticResult');xlabel('f');ylabel('A(f)');subplot(2,1,2);plot(ft,180/pi*ph,'r-',ft,180/pi*angle(B1),'b-.');title('PhaseofSweepSignal');legend('NumericResult','TheoreticResult');xlabel('f');ylabel('\Psi(f)');------------------------------------------------------------------functiony=fsin(x,ot,yd)a=x(1);%amplitudeotisomega*tph=x(2);%phaseinradiany=a*sin(ot+ph)-yd;-------------------------------------------------------------------由于相隔此次的频率相距很近，因而把上一次拟合的结果作为本次的初值，不但可以保证初始值给得非常合理，同时可以加快计算速度。另外要强调的是尽管如此，由于每个频率段都要使用拟合，因而分段曲线拟合方法计算速度比较慢。正弦扫频信号幅值及相位的提取(4)方法3分段两点求解在[f,f+df]区间内，利用两点求出两个未知数A，ψ，在[f,f+df]区间内对A，ψ取平均。图5给出了有相位差时曲线拟合提取的幅值与相位同理论结果的比较，从图中可以看出计算结果与真实值基本重合。由于算例中的扫频信号是理想的正弦扫频信号，因而两点求解能够精确计算得到真实值。8图5使用分段两点求解提取的频域幅值、相位下面给出了实现分段两点求解提取扫频信号的频域幅值、相位的Matlab代码。---------------------------------------------------------------------%Decomposetheamplitudeandphasefromthesweepsignal%Gettingtwoparametersthroughtwopointsf1=5;%theinitialfreqs=4;%sweepratefr=50;%Resonantfreqaf=[];%amplitudephf=[];%phasek1=0.02;%dampingratiodf=0.01;%freqincreaseforfa=40:df:60t1=60/s*log2(fa/f1);9t2=60/s*log2((fa+df)/f1);ta=t1:0.001:t2;N=length(ta);ft=f1*2.^(s/60*ta);A1=sin(2*pi*ft.*ta);lamb=ft/fr;B1=1./(1-lamb.^2+j*2*k1*lamb);%transferfunctionA2=abs(B1).*sin(2*pi*ft.*ta+angle(B1));amp=0;ph=0;fornk=2:NT1=[sin(2*pi*ft(nk-1)*ta(nk-1)),cos(2*pi*ft(nk-1)*ta(nk-1));...sin(2*pi*ft(nk)*ta(nk)),cos(2*pi*ft(nk)*ta(nk))];b1=[A2(nk-1);A2(nk)];X=T1\b1;amp=amp+sqrt(sum(X.^2));ph=ph+angle(X(1)+j*X(2));endaf=[af,amp/(N-1)];phf=[phf,ph/(N-1)];endfa=40:df:60;lamb=fa/fr;bf=abs(1./(1-lamb.^2+j*2*k1*lamb));subplot(2,1,1);plot(fa,af,'r-',fa,bf,'b-.');legend('NumericResult','TheoreticResult');title('MagnitudeoftheSweepSignal');xlabel('f');ylabel('A(f)');10subplot(2,1,2);plot(fa,180/pi*phf,'r-',fa,180/pi*angle(1./(1-lamb.^2+j*2*k1*lamb)),'b-.');legend('NumericResult','TheoreticResult');title('PhaseoftheSweepSignal');xlabel('f');ylabel('\Psi(f)');--------------------------------------------------------------------尽管分段两点求解计算精度高，求解速度快，但是由于在计算两参数时使用了矩阵求逆(2×2的矩阵)，因而可能会由于相邻两点的线性相关导致矩阵退化计算无法进行而中止。正弦扫频信号幅值及相位的提取(5)方法4峰值包络