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

正弦扫频信号幅值及相位的提取(1)正弦振动控制系统提供输入的扫频信号 ，对于对数扫频， ，其中Sr为对数扫描率，若频响函数为 则系统输出为 。 测量系统中可得到Calo信号及响应信号，通过对二者进行数据处理，可得到频域下的响应 。不知道LMS的信号采集软件是如何提取频域响应的，个人认为软件计算速度有限，LMS应该是通过硬件实现的。下面我提供几种方法并进行比较。 算例对于Calo信号 ，频响函数为 ，其中 ，信号采样率为1000次/秒，图1给出了时域下的响应信号。图 1 时域下的响应信号正弦扫频信号幅值及相位的提取(2)方法1 分段FFT 在f, f+df区间内对Calo信号、响应信号进行FFT变换，二者在频率f处的谱值比即为频响函数在f处的值。此方法的缺陷是由于信号采样率为1000Hz，而f, f+df的区间很窄，在此区间下时域的点不会很多，因而FFT的频率分辨率不高。 对于没有相位差的扫频信号，此方法能较好的提取幅值。图2给出了使用此方法提取的幅值与理论结果比较，由图中可以看出二者基本吻合。图 2 使用分段FFT提取的频域幅值 对于有相位差的扫频信号，则要对结果进行光滑处理，Matlab的smooth函数提供了这一功能。图3给出了有相位差时分段FFT提取的幅值与相位同理论结果的比较，从图中可以看出在频域峰值处分段FFT比理论值大，在其余频段二者吻合较好。图 3 使用分段FFT提取的频域幅值、相位 下面给出了实现分段FFT提取扫频信号的频域幅值、相位的Matlab代码。-% Decompose the amplitude and phase from the sweep signal% Local fft and smooth are employed.f1=5; % the initial freqs=4; % sweep ratefr=50; % Resonant freqaf=; % amplitudepf=; % phasek1=0.02; % damping ratiodf=0.01; % freq intervalfor fa=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); A1=sin(2*pi*ft.*ta); lamb=ft/fr; B1=1./(1-lamb.2+j*2*k1*lamb); % transfer function A2=abs(B1).*sin(2*pi*ft.*ta+angle(B1); ffreq=exp(-j*2*pi*(fa-400)*ta); % freq shift for time domain spa=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); % Key trickfa=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(Numeric Result,Theoretic Result);title(Amplitude of Sweep Signal);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(Numeric Result,Theoretic Result);title(Phase of Sweep Signal);xlabel(f);ylabel(Psi(f);- 分段FFT提取方法计算速度一般，不会出现异常而中止，计算精度基本也能保证。正弦扫频信号幅值及相位的提取(3)方法2 分段曲线拟合 在f, f+df区间内，假定A，不变，此区间内在时域内对其拟合。图4给出了有相位差时曲线拟合提取的幅值与相位同理论结果的比较，从图中可以看出计算结果与真实值吻合非常好。图 4 使用分段曲线拟合提取的频域幅值、相位 分段曲线拟合提取的结果精度非常高，但是由于是拟合方法，因而可能会由于初始值给的不合理或拟合关系式不恰当而出现迭代次数超过规定值从而导致计算中止。下面给出了实现分段曲线拟合提取扫频信号的频域幅值、相位的Matlab代码。-% Decompose the amplitude and phase from the sweep signal% Local curve fit is applied.f1=5; % the initial freqs=4; % sweep ratefr=50; % Resonant freqaf=; % amplitudek1=0.02; % damping ratiodf=0.01; % freq increasea=; % amplitudeph=; % phasex0=1,0; % initial guessfor fa=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); % transfer function A1=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; % Key trick a=a,x(1); ph=ph,x(2);endft=40:df:60;lamb=ft/fr;B1=1./(1-lamb.2+j*2*k1*lamb); % transfer functionsubplot(2,1,1);plot(ft,a,r-,ft,abs(B1),b-.);title(Magnitude of Sweep Signal);legend(Numeric Result,Theoretic Result);xlabel(f);ylabel(A(f);subplot(2,1,2);plot(ft,180/pi*ph,r-,ft,180/pi*angle(B1),b-.);title(Phase of Sweep Signal);legend(Numeric Result,Theoretic Result);xlabel(f);ylabel(Psi(f);-function y=fsin(x,ot,yd)a=x(1); % amplitude ot is omega*tph=x(2); % phase in radiany=a*sin(ot+ph)-yd;- 由于相隔此次的频率相距很近，因而把上一次拟合的结果作为本次的初值，不但可以保证初始值给得非常合理，同时可以加快计算速度。另外要强调的是尽管如此，由于每个频率段都要使用拟合，因而分段曲线拟合方法计算速度比较慢。正弦扫频信号幅值及相位的提取(4)方法3 分段两点求解 在f, f+df区间内，利用两点求出两个未知数A，在f, f+df区间内对A，取平均。图5给出了有相位差时曲线拟合提取的幅值与相位同理论结果的比较，从图中可以看出计算结果与真实值基本重合。由于算例中的扫频信号是理想的正弦扫频信号，因而两点求解能够精确计算得到真实值。图 5 使用分段两点求解提取的频域幅值、相位 下面给出了实现分段两点求解提取扫频信号的频域幅值、相位的Matlab代码。-% Decompose the amplitude and phase from the sweep signal% Getting two parameters through two pointsf1=5; % the initial freqs=4; % sweep ratefr=50; % Resonant freqaf=; % amplitudephf=; % phasek1=0.02; % damping ratiodf=0.01; % freq increasefor fa=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); A1=sin(2*pi*ft.*ta); lamb=ft/fr; B1=1./(1-lamb.2+j*2*k1*lamb); % transfer function A2=abs(B1).*sin(2*pi*ft.*ta+angle(B1); amp=0; ph=0; for nk=2:N T1=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=T1b1; amp=amp+sqrt(sum(X.2); ph=ph+angle(X(1)+j*X(2); end af=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(Numeric Result,Theoretic Result);title(Magnitude of the Sweep Signal);xlabel(f);ylabel(A(f);subplot(2,1,2);plot(fa,180/pi*phf,r-,fa,180/pi*angle(1./(1-lamb.2+j*2*k1*lamb),b-.);legend(Numeric Result,Theoretic Result);title(Phase of the Sweep Signal);xlabel(f);ylabel(Psi(f);- 尽管分段两点求解计算精度高，求解速度快，但是由于在计算两参数时使用了矩阵求逆(22的矩阵)，因而可能会由于相邻两点的线性相关导致矩阵退化计算无法进行而中止。正弦扫频信号幅值及相位的提取(5)方法4 峰值包络 提取时域扫频曲线的峰值与谷值，通过对Colo信号、响应信号进行分析得到相应的频率及相位信息。图6给出了峰值包络提取的幅值与相位同理论结果的比较。图中进行了光滑处理，从图中可以看出提取的幅值在各频率都低于精确值，由于振幅是随频率变化的，时域的峰值只是振幅与相位二者综合后形成的，即该处出现峰值相位并非为/2，峰值只是幅值与相位的正弦之积，因而实际幅值比峰值大。由于相位是间接提取的，相比于幅值误差更大。图 6 使用峰值包络提取的频域幅值、相位 下面给出了利用峰值包络提取扫频信号的频域幅值、相位的Matlab代码。这里可能会出现当精确值为-180由于误差导致小于-180，由于相位一般定义为-180,180，因而此时相位接近180，尽管只是表示形式的差别，但显示在相位曲线上就谬以千里了，故这里进行了特殊处理，即将相位表示为精确值与误差值之和，由于相位的误差值在很小的范围内，因而就避免了上述相位表示的问题。图中当相位低于-180时，即以小于-180的形式表示。-% Decompose the amplitude and phase from the sweep signal% Searching peaks from time domain curvef1=5; % the initial freqs=4; % sweep ratefr=50; % Resonant freqaf=; % amplitudepf=; % phasek1=0.02; % damping ratiodf=20; % freq increasenk=; % the number of peak points in the curvefa=40;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);A1=sin(2*pi*ft.*ta);lamb=ft/fr;B1=1./(1-lamb.2+j*2*k1*lamb); % transfer functionA2=abs(B1).*sin(2*pi*ft.*ta+angle(B1);for i=2:N-1 if A2(i)A2(i-1)&A2(i)A2(i+1) % peak nk=nk,i; pa=angle(sin(2*pi*ft(i)*ta(i)+j*cos(2*pi*ft(i)*ta(i); pf=pf,pa; elseif A2(i)A2(i-1)&A2(i)A2(i+1) % valley nk=nk,i; pa=angle(-sin(2*pi*ft(i)*ta(i)-j*cos(2*pi*ft(i)*ta(i); pf=pf,pa; endendtb=ta(nk);fa=f1*2.(s/60*tb); % the corresponding freq for the peakaf=abs(A2(nk);frb=fa/fr;pf2=angle(1./(1-frb.2+j*2*k1*frb);pf=pf2+angle(cos(pf-pf2)+j*sin(pf-pf2); %