9种谱校正方法及matlab代码
9种谱校正方法
9种谱校正方法及matlab 程序代码采样间隔归一化成1T ∆=,采样长度为N .这样FFT 离散谱线为0,1)i X i N =-（,相应的频率分辨率2/(1/)N f N ωπ∆=∆=. 设FFT 离散谱线局部极高谱线为m (为了数学上简洁,假定从0开始,注意在MA TLAB 环境下数组实际操作的是从1开始),记频偏量δωδω=∆. 我们需要使用谱线m 和与之相邻一条次高谱线，记这连续两条谱线中左边一条序号为M (当次高谱线在m 左侧时1M m =-,反之M m =).下面列出若干算法的δ计算公式1. 加矩形窗的精确谱校正[1]i i i X U jV =+111()sin()()cos()M M M M opt M MV V M U U M K U U ωω+++-∆+-∆=- 1211cos()sin()cos()sin()opt M M opt M M K M Z V U M K M Z V U M ωωωωωω++-∆⎡⎤=+⎢⎥∆⎣⎦-∆+∆⎡⎤=+⎢⎥∆+∆⎣⎦ 2121cos()cos()()Z M Z M M m Z Z ωωωδ∆+∆-∆=+-- 2. 加矩形窗情形,采用解析单频模型的幅值比校正[1, 2]11||()||||M M M X M m X X δ++=+-+ 3. 加汉宁窗情形,采用解析单频模型的幅值比校正[1, 2]112||||()||||M M M M X X M m X X δ++-=+-+ 4. 加矩形窗情形,采用解析单频模型的复比值校正[3]11Re ()M M M X M m X X δ++⎛⎫=+- ⎪-⎝⎭5. 加汉宁窗情形,采用解析单频模型的复比值校正[3]112()M M M MX X M m X X δ+++=+-- 6. 加矩形窗情形,采用解析单频模型的复合复比值校正[3]11Re ()M m M M X M m X X δ++⎛⎫=+- ⎪-⎝⎭ 11m R m mX X X δ++=-,1111m m L m m m m X X X X X X δ---=-=-- 0.5)0.5)m L m R δδδδδ=-++（（7. 加汉宁窗情形,采用解析单频模型的复合复比值校正[3]112Re ()M M m M M X X M m X X δ++⎛⎫+=+-⎪-⎝⎭ 112m m R m mX X X X δ+++=-,1111221m m m m L m m m m X X X X X X X X δ----++=-=-- 0.5)0.5)m L m R δδδδδ=-++（（8. 加矩形窗,Quin 校正[4]11Re()Re(),Re()Re()m m L R m m X X X X αα-+== 11L R L R L Rααδδαα==---，, 00 R R L R δδδδδ>>⎧=⎨⎩当且其它9. 加汉宁窗,Quin 校正[4]11Re()Re(),Re()Re()m m L R m m X X X X αα-+== 212111L RL R L Rααδδαα++==---，,00 R R L R δδδδδ>>⎧=⎨⎩当且其它References1. Schoukens, J., R. Pintelon,H. Van Hamme. The interpolated fast Fourier transform: Acomparative study . IEEE Transactions on Instrumentation and Measurement. 1992, 41(2):226-232.2. 谢明,丁康. 频谱分析的校正方法.振动工程学报. 1994, 7(2): 172-179.3. 陈奎孚, 王建立,张森文. 频谱校正的复比值法.振动工程学报(已投). 2007.4. Quinn, B.G. Estimating frequency by interpolation using Fourier coefficients.IEEETransactions on Signal Processing. 1994, 42(5): 1264-1268.%========================这是调用调试==================DT=1;N=1024;PHI=pi/3;Ampl=1;CiR=11.9; %cycles in recordFreq=CiR/(DT*N); %frequencyTV=[0:N-1];DatVec=Ampl*cos(Freq*TV*2*pi+PHI);FV=fft(DatVec);figuresubplot(2,1,1);plot(TV,DatV ec);subplot(2,1,2);plot(abs(FV(1:round(N/2.56))));grid on[MV,MI]=max(abs(FV));%加矩形窗的解析校正--1FreqShift=SpecCorr(FV,MI,N,1);%加矩形窗的解析单频模型校正--2FreqShift=SpecCorr(FV,MI,N,2);%加汉宁窗的解析单频模型校正--3HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,3);%加矩形窗的解析单频模型校正+复比值法--4FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,4);%加汉宁窗的解析单频模型校正+复比值法--5HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,5);%加矩形窗的解析单频模型校正+复比值法+左右平均--6FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,6);%加汉宁窗的解析单频模型校正+复比值法+左右平均--7HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,7);%加矩形窗的解析单频模型校正+Quinn算法--8FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,8);%加汉宁窗的解析单频模型校正+Quinn算法--9HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,9);===========这是子程序======================%spectrum correction assemble% the sampling interval is 1 s (or unitary)%Input: SpecVec--Discrte Fourier Spectrum from FFT% PeakIdx--the peak index, noting the matrix in MatLab start from 1% TL--the length (or the point number) of the FFT% method--correction method% output: PeakShift-- the corrected peak shifting from the peak in discrete% spectrumfunction PeakShift=SpecCorr(SpecVec,PeakIdx,TL,method)% picking up the second highest spectrum lineif(abs(SpecVec(PeakIdx-1))>abs(SpecVec(PeakIdx+1)))IP=[PeakIdx-1,PeakIdx];ShiftCorr=-1; %shift aligning with the PeakIdxelseIP=[PeakIdx,PeakIdx+1];ShiftCorr=0; %shift aligning with the PeakIdxendII=IP(1)-1; % noting that the index of a matrix in MATLAB starts from 1, not zero if(method==1) %an analyitic solution for rectangular windowU=real(SpecVec(IP));V=imag(SpecVec(IP));DW=2*pi/TL;KOPT=(sin(II*DW)*(V(2)-V(1))+cos(II*DW)*(U(2)-U(1)))/(U(2)-U(1));Z=V.*(KOPT-cos((IP-1)*DW))./(sin(DW*(IP-1)))+U;Tmp1=(Z(2)*cos(DW*(II+1))-Z(1)*cos(DW*II))/(Z(2)-Z(1));PeakPos=acos(Tmp1)/DW;PeakShift=PeakPos-(PeakIdx-1);elseif(method==2) %based on the analytical-single-tone model for rectangular window PeakShift=abs(SpecVec(IP(2)))/(abs(SpecVec(IP(2)))+abs(SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==3) %based on the analytical-single-tone model for Hanning windowPeakShift=(2*abs(SpecVec(IP(2)))-abs(SpecV ec(IP(1))))/(abs(SpecV ec(IP(2)))+abs(SpecVec(IP(1 ))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==4) %based on the analytical-single-tone model for rectangular window with complex correctionPeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==5) %based on the analytical-single-tone model for Hanning window with complex correctionPeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1)));PeakShift=real(PeakShift)+ShiftCorr; %shift aligning with the PeakIdx elseif(method==6) %based on the analytical-single-tone model for rectangular window with complex correction+averagePeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));MaxPeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxLeftShift=real(SpecVec(PeakIdx)/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1)))-1;RightShift=real(SpecVec(PeakIdx+1)/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx)));%averagePeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift;elseif(method==7) %based on the analytical-single-tone model for Hanning window with complex correction+average,????PeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1)));MaxPeakShift=real(PeakShift)+ShiftCorr; %shift aligning with the PeakIdxLeftShift=(2*SpecVec(PeakIdx)+SpecV ec(PeakIdx-1))/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1))-1;RightShift=(2*SpecVec(PeakIdx+1)+SpecV ec(PeakIdx))/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx ));%averagePeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift;elseif(method==8) %Quinn method for the rectangular windowa1 = real(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %lefta2 = real(SpecV ec(PeakIdx+1)/SpecVec(PeakIdx)); %rightLeftShift = a1/(1-a1);RightShift = -a2/(1-a2);if (LeftShift>0 & RightShift>0)PeakShift = RightShift;elsePeakShift = LeftShift;endelseif(method==9) %Quinn method for the Hanning windowa1 = real(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %lefta2 = real(SpecV ec(PeakIdx+1)/SpecVec(PeakIdx)); %rightLeftShift = (2*a1+1)/(1-a1);RightShift = -(2*a2+1)/(1-a2);if (LeftShift>0 & RightShift>0)PeakShift = RightShift;elsePeakShift = LeftShift;endendend。
谱线弯曲校正的方法
校正谱线弯曲的方法1、利用棱镜-光栅（P-G）组合分光元件并结合系统物镜畸变校正谱线弯曲.原理：分别计算了棱镜和光栅产生的谱线弯曲以及P-G组合元件产生的光谱弯曲，分析了棱镜和光栅的谱线弯曲特性，并基于此设计了P-G组合分光元件和消谱线弯曲成像光谱仪结构。
所设计的成像光谱仪谱线弯曲和光谱弯曲均小于2微米。
光谱分辨率在2纳米以上，系统像散较小，成像质量满足空间分辨率要求。
设计结果表明，采用P-G组合分光元件校正中心波长谱线弯曲，使剩余谱线弯曲和光谱弯曲对称分布，然后结合准直物镜和成像物镜的畸变补偿校正剩余的谱线弯曲和光谱弯曲，可实现全工作谱段的谱线弯曲校正。
方法:引入棱镜和光栅谱线弯曲公式，分析了棱镜和光栅谱线弯曲的相关性质，并推导了P-G组合分光元件校正谱线弯曲的公式，由此建立校正谱线弯曲的方案。
1、 棱镜谱线弯曲图1为棱镜折射产生谱线弯曲的示意图。
MQ为狭缝半高，Q为狭缝中心，QO为系统光轴。
狭缝中心Q点发出的光线经过棱镜的主截面ABC;狭缝的上端点M发出的光线MO通过棱镜的截面ABD,由于两个截而的色散顶角不相等，使得两点经过棱镜的色散角不同，成像后谱线发生弯曲。
假设系统中准直物镜和成像物镜焦距相同，经棱镜色散后的谱线弯曲公式可表示2、 光栅谱线弯曲光栅谱线弯曲的原理与棱镜相似，狭缝上不同点发出的光线入射在光栅时相对于主截而是倾斜的，其衍射光线在主平而的投影也将不同，造成衍射光成像时谱线弯曲。
假设系统中准直物镜和成像物镜焦距相同，透射平而光栅的谱线弯曲公式可表示为比较两式，可知棱镜和光栅的谱线弯曲方向相反，在设计过程中可利用棱镜和光栅补偿校正谱线弯曲。
3、 校正中心波长谱线弯曲计算图为中心波长光通过分光元件的光路示意图，狭缝中心发出的光入射到棱镜的表面#经棱镜折射后进入体全息相位光栅。
结合折射定律及上述两种谱线弯曲公式可得出求解上式组成的方程组可求得棱镜两面的倾斜角从而确定P-G组合分光元件的结构。
MATLAB技术图像矫正方法
MATLAB技术图像矫正方法MATLAB技术在图像矫正中的应用引言：图像是我们日常生活中非常普遍的一种媒介，在各个领域都得到了广泛的应用。
然而，由于摄像设备的限制以及环境因素的影响，图像中往往会存在各种各样的畸变现象。
为了获得更加真实准确的图像，我们需要对这些畸变进行矫正。
本文将介绍MATLAB技术中常用的图像矫正方法。
一、几何矫正方法几何矫正方法主要用于纠正图像中由于成像设备本身造成的畸变，如镜头畸变、透视畸变等。
其中，常用的方法有摄像机标定和透视变换。
1. 摄像机标定摄像机标定是一种通过获取摄像机参数来矫正图像畸变的方法。
它通过拍摄特定的棋盘格图案，在摄像机视野范围内进行多个角度和位置的变换，从而计算出摄像机的内参和外参。
通过这些参数，可以将图像进行几何矫正，消除镜头畸变。
MATLAB提供了内置的相机标定工具箱，可以方便地进行摄像机的标定。
它提供了棋盘格图案的自动检测、校正板的自动提取等功能，大大简化了标定的过程。
2. 透视变换透视变换是一种通过投影变换矫正图像透视畸变的方法。
透视畸变是由于成像平面和实际场景之间的角度关系引起的，常见于拍摄建筑物或者远处景物时。
透视变换可以将图像中的平行线在成像平面上呈现为平行线，从而消除透视畸变。
在MATLAB中，可以利用仿射变换函数和透视变换函数实现透视矫正。
通过选取图像中的多个关键点，计算出变换矩阵，再进行透视变换，即可实现图像的透视矫正。
二、亮度矫正方法亮度矫正方法用于调整图像亮度和对比度，使图像更加清晰明亮。
其中，常用的方法有直方图均衡化和灰度拉伸。
1. 直方图均衡化直方图均衡化是一种通过增强图像对比度来提高图像亮度的方法。
它通过分布图像中像素值的频率来调整图像灰度级的分布情况，使其更加均匀。
这样可以使得图像的细节更加清晰，颜色更加鲜艳。
MATLAB中提供了直方图均衡化的函数，只需输入待矫正的图像，就可以得到均衡化后的图像。
同时，还可以通过调整直方图均衡化函数的参数，进一步调整亮度和对比度的效果。
matlab几何纠正，间接法，双线性内插
matlab几 何 纠 正 ， 间 接 法 ， 双 线 性 内 插
超简洁，超级快，两个文件 datapre.m文件代码: global X; global Y; global A; global l; global i; global I; global m; global n; global k; global dx;
['待纠正图像Y: ',num2str(Y)]}; i=i+1; if(rem(i,2)==1) k=(i+1)/2; A(k,1)=1; A(k,2)=X; A(k,3)=Y; else
l(i/2,1)=X; l(i/2,2)=Y; end [m,n]=size(l); if(m==4) datacursormode off; dx=inv(A'*A)*(A'*l); [m,n,k]=size(I); dis=sqrt((A(2,2)-A(1,2))*(A(2,2)-A(1,2))+(A(2,3)-A(1,3))*(A(2,3)-A(1,3))); dist=sqrt((l(2,1)-l(1,1))*(l(2,1)-l(1,1))+(l(2,2)-l(1,2))*(l(2,2)-l(1,2))); sca=dist/dis; cosa=(dx(2,1)+dx(3,2))/(2*sca); sina=(-dx(2,2)+dx(3,1))/(2*sca); alfa=atan(sina/cosa); alfa=alfa*180/3.1415926; ResImage=imresize(I,sca,'bilinear'); ResImage=imrotate(ResImage,alfa,'bilinear');
光谱校正的校正方法
光谱校正的校正方法
光谱校正是一种用于纠正光谱数据中由于仪器或实验条件引起的偏差的方法。
以下是一些常见的光谱校正方法：
1. 波长校正：通过测量已知标准物质的光谱特征峰位置，对光谱数据进行波长校正，以确保光谱的波长准确性。
2. 强度校正：通过测量已知标准物质的光谱强度，对光谱数据进行强度校正，以确保光谱的强度准确性。
3. 基线校正：通过减去光谱数据中的背景信号或基线，对光谱数据进行基线校正，以去除背景干扰。
4. 光谱平滑：通过对光谱数据进行平滑处理，减少噪声和光谱波动，以提高光谱的质量。
5. 光谱归一化：通过对光谱数据进行归一化处理，使得不同光谱之间具有可比性。
6. 光谱去卷积：通过对光谱数据进行去卷积处理，去除光谱中的卷积效应，以提高光谱的分辨率。
这些光谱校正方法可以单独或组合使用，以满足不同的应用需求。
在进行光谱校准时，需要根据具体情况选择合适的校正方法，并进行适当的参数设置和优化，以确保校正的准确性和有效性。
求信号功率谱时候用下面不同方法
求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？功率谱密度的幅植的具体意义是什么？？下面是一些不同方法计算同一信号的matlab 程序！欢迎大家给点建议！直接法：直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
Matlab代码示例：clear;Fs=1000; %采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));window=boxcar(length(xn)); %矩形窗nfft=1024;[Pxx,f]=periodogram(xn,window,nfft,Fs); %直接法plot(f,10*log10(Pxx));间接法：间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
Matlab代码示例：clear;Fs=1000; %采样频率n=0:1/Fs:1;%产生含有噪声的序列xn=cos(2*pi*40*n)+3*cos(2*pi*100*n)+randn(size(n));nfft=1024;cxn=xcorr(xn,'unbiased'); %计算序列的自相关函数CXk=fft(cxn,nfft);Pxx=abs(CXk);index=0:round(nfft/2-1);k=index*Fs/nfft;plot_Pxx=10*log10(Pxx(index+1));plot(k,plot_Pxx);改进的直接法：对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
matlab中的谱方法
在MATLAB中，谱方法通常用于信号处理、频谱分析、滤波以及其他与频域相关的操作。
以下是一些常见的MATLAB函数和工具，用于实现谱方法：1. **傅立叶变换**：MATLAB提供了`fft`函数，用于计算信号的快速傅立叶变换（FFT）。
它允许你将信号从时域转换到频域。
```matlabX = fft(x);```2. **功率谱密度**：使用谱方法来估计信号的功率谱密度（PSD）。
`pwelch`和`periodogram`是两个常用的函数，用于估计信号的功率谱密度。
```matlab[Pxx, f] = pwelch(x, window, overlap, nfft, fs);```3. **滤波**：使用谱方法来设计和应用数字滤波器，以对信号进行滤波。
MATLAB中有一些滤波函数，如`filter`和`designfilt`。
```matlaby = filter(b, a, x);```4. **频域可视化**：使用`plot`等函数可以可视化频域数据，以便分析信号的频谱内容。
```matlabplot(f, 10*log10(Pxx));xlabel('Frequency (Hz)');ylabel('Power/Frequency (dB/Hz)');```5. **信号合成**：你可以使用逆傅立叶变换将频域信号合成回时域信号。
```matlabx_reconstructed = ifft(X);```这些是MATLAB中常见的一些谱方法的示例。
你可以根据你的具体需求和信号处理任务来选择合适的工具和函数。
MATLAB的文档和示例也可以提供更多帮助和指导。
自动控制原理课件控制系统校正MATLAB程序
表6-1
相位滞后校正
常见无源校正装置
相位超前校正 相位滞后-超前校正
RC
网络
R1 R1 R2 1 R1C1 1 R1C1 (1s 1)( 2 s 1) K (1s 1) K R1 2s 1 R1 G(s) R R R G ( s ) 1 2 1K G1 (s) R2C2 (1s 1)( 2 s 1) 2 C1R1 R2 2R 1C2 s K 2s 1 1s 1 R1 R2 R1R R R 1 22 R C 1 1 1 (1s 1)( 2 s 1) 1 ( R1 R2 )C2 1 2 R C R C 1 1 1 ≥ 1 2 1 1 1 R1R2 ( s 1)( s 1) ( R R ) C 1 R 2 2 C1 1 1R1C1 2 2 R2C2 1 1 ( R1 2 2 2 )C R R R R R R 2 2 2 1 C1 1 1 2 2 C1 2 C2 2 RR 2 2 R2C2 2 1 2 C2 1 R1C1 1 ≥ 2 R1 R 2 R 1 R2 2 2 1 1 1 2 1≥ 2 1≥ 2 2 R2C2 K
(2) 校正后系统性能分析：
设校正装置的传递函数为 Gc(s) = K(τds＋1)，为了说明 相位超前校正对系统性能的影响，取 K = 1 微分时间常数取 τd=T1 = 0.2s(抵消大惯性环节的相位滞后影响)，则系统的开 环传递函数变为：
K K1 G G o( s ) G1 ( s )Gc ( s ) K c ( d s 1) s (T1s 1)(T2 s 1) s (T2
6.2.3
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
9种谱校正方法及matlab 程序代码采样间隔归一化成1T ∆=,采样长度为N .这样FFT 离散谱线为0,1)i X i N =-（,相应的频率分辨率2/(1/)N f N ωπ∆=∆=. 设FFT 离散谱线局部极高谱线为m (为了数学上简洁,假定从0开始,注意在MATLAB 环境下数组实际操作的是从1开始),记频偏量δωδω=∆. 我们需要使用谱线m 和与之相邻一条次高谱线，记这连续两条谱线中左边一条序号为M (当次高谱线在m 左侧时1M m =-,反之M m =).下面列出若干算法的δ计算公式1. 加矩形窗的精确谱校正[1]i i i X U jV =+111()sin()()cos()M M M M opt M MV V M U U M K U U ωω+++-∆+-∆=- 1211cos()sin()cos()sin()opt M M opt M M K M Z V U M K M Z V U M ωωωωωω++-∆⎡⎤=+⎢⎥∆⎣⎦-∆+∆⎡⎤=+⎢⎥∆+∆⎣⎦ 2121cos()cos()()Z M Z M M m Z Z ωωωδ∆+∆-∆=+-- 2. 加矩形窗情形,采用解析单频模型的幅值比校正[1, 2]11||()||||M M M X M m X X δ++=+-+ 3. 加汉宁窗情形,采用解析单频模型的幅值比校正[1, 2]112||||()||||M M M M X X M m X X δ++-=+-+ 4. 加矩形窗情形,采用解析单频模型的复比值校正[3]11Re ()M M M X M m X X δ++⎛⎫=+- ⎪-⎝⎭5. 加汉宁窗情形,采用解析单频模型的复比值校正[3]112()M M M MX X M m X X δ+++=+-- 6. 加矩形窗情形,采用解析单频模型的复合复比值校正[3]11Re ()M m M M X M m X X δ++⎛⎫=+- ⎪-⎝⎭ 11m R m mX X X δ++=-,1111m m L m m m m X X X X X X δ---=-=-- 0.5)0.5)m L m R δδδδδ=-++（（7. 加汉宁窗情形,采用解析单频模型的复合复比值校正[3]112Re ()M M m M M X X M m X X δ++⎛⎫+=+-⎪-⎝⎭ 112m m R m mX X X X δ+++=-,1111221m m m m L m m m m X X X X X X X X δ----++=-=-- 0.5)0.5)m L m R δδδδδ=-++（（8. 加矩形窗,Quin 校正[4]11Re()Re(),Re()Re()m m L R m m X X X X αα-+== 11L R L R L Rααδδαα==---，, 00 R R L R δδδδδ>>⎧=⎨⎩当且其它9. 加汉宁窗,Quin 校正[4]11Re()Re(),Re()Re()m m L R m m X X X X αα-+== 212111L RL R L Rααδδαα++==---，,00 R R L R δδδδδ>>⎧=⎨⎩当且其它References1. Schoukens, J., R. Pintelon,H. Van Hamme. The interpolated fast Fourier transform: Acomparative study . IEEE Transactions on Instrumentation and Measurement. 1992, 41(2):226-232.2. 谢明,丁康. 频谱分析的校正方法.振动工程学报. 1994, 7(2): 172-179.3. 陈奎孚, 王建立,张森文. 频谱校正的复比值法.振动工程学报(已投). 2007.4. Quinn, B.G. Estimating frequency by interpolation using Fourier coefficients.IEEETransactions on Signal Processing. 1994, 42(5): 1264-1268.%========================这是调用调试==================DT=1;N=1024;PHI=pi/3;Ampl=1;CiR=11.9; %cycles in recordFreq=CiR/(DT*N); %frequencyTV=[0:N-1];DatVec=Ampl*cos(Freq*TV*2*pi+PHI);FV=fft(DatVec);figuresubplot(2,1,1);plot(TV,DatV ec);subplot(2,1,2);plot(abs(FV(1:round(N/2.56))));grid on[MV,MI]=max(abs(FV));%加矩形窗的解析校正--1FreqShift=SpecCorr(FV,MI,N,1);%加矩形窗的解析单频模型校正--2FreqShift=SpecCorr(FV,MI,N,2);%加汉宁窗的解析单频模型校正--3HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,3);%加矩形窗的解析单频模型校正+复比值法--4FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,4);%加汉宁窗的解析单频模型校正+复比值法--5HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,5);%加矩形窗的解析单频模型校正+复比值法+左右平均--6FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,6);%加汉宁窗的解析单频模型校正+复比值法+左右平均--7HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,7);%加矩形窗的解析单频模型校正+Quinn算法--8FV=fft(DatVec);FreqShift=SpecCorr(FV,MI,N,8);%加汉宁窗的解析单频模型校正+Quinn算法--9HanDat=DatVec.*hanning(N,'periodic')';FV=fft(HanDat);FreqShift=SpecCorr(FV,MI,N,9);===========这是子程序======================%spectrum correction assemble% the sampling interval is 1 s (or unitary)%Input: SpecVec--Discrte Fourier Spectrum from FFT% PeakIdx--the peak index, noting the matrix in MatLab start from 1% TL--the length (or the point number) of the FFT% method--correction method% output: PeakShift-- the corrected peak shifting from the peak in discrete% spectrumfunction PeakShift=SpecCorr(SpecVec,PeakIdx,TL,method)% picking up the second highest spectrum lineif(abs(SpecVec(PeakIdx-1))>abs(SpecVec(PeakIdx+1)))IP=[PeakIdx-1,PeakIdx];ShiftCorr=-1; %shift aligning with the PeakIdxelseIP=[PeakIdx,PeakIdx+1];ShiftCorr=0; %shift aligning with the PeakIdxendII=IP(1)-1; % noting that the index of a matrix in MATLAB starts from 1, not zero if(method==1) %an analyitic solution for rectangular windowU=real(SpecVec(IP));V=imag(SpecVec(IP));DW=2*pi/TL;KOPT=(sin(II*DW)*(V(2)-V(1))+cos(II*DW)*(U(2)-U(1)))/(U(2)-U(1));Z=V.*(KOPT-cos((IP-1)*DW))./(sin(DW*(IP-1)))+U;Tmp1=(Z(2)*cos(DW*(II+1))-Z(1)*cos(DW*II))/(Z(2)-Z(1));PeakPos=acos(Tmp1)/DW;PeakShift=PeakPos-(PeakIdx-1);elseif(method==2) %based on the analytical-single-tone model for rectangular window PeakShift=abs(SpecVec(IP(2)))/(abs(SpecVec(IP(2)))+abs(SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==3) %based on the analytical-single-tone model for Hanning windowPeakShift=(2*abs(SpecVec(IP(2)))-abs(SpecV ec(IP(1))))/(abs(SpecV ec(IP(2)))+abs(SpecVec(IP(1 ))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==4) %based on the analytical-single-tone model for rectangular window with complex correctionPeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));PeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxelseif(method==5) %based on the analytical-single-tone model for Hanning window with complex correctionPeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1)));PeakShift=real(PeakShift)+ShiftCorr; %shift aligning with the PeakIdx elseif(method==6) %based on the analytical-single-tone model for rectangular window with complex correction+averagePeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));MaxPeakShift=PeakShift+ShiftCorr; %shift aligning with the PeakIdxLeftShift=real(SpecVec(PeakIdx)/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1)))-1;RightShift=real(SpecVec(PeakIdx+1)/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx)));%averagePeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift;elseif(method==7) %based on the analytical-single-tone model for Hanning window with complex correction+average,????PeakShift=(2*SpecVec(IP(2))+SpecVec(IP(1)))/(SpecVec(IP(2))-SpecVec(IP(1)));MaxPeakShift=real(PeakShift)+ShiftCorr; %shift aligning with the PeakIdxLeftShift=(2*SpecVec(PeakIdx)+SpecV ec(PeakIdx-1))/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1))-1;RightShift=(2*SpecVec(PeakIdx+1)+SpecV ec(PeakIdx))/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx ));%averagePeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift;elseif(method==8) %Quinn method for the rectangular windowa1 = real(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %lefta2 = real(SpecV ec(PeakIdx+1)/SpecVec(PeakIdx)); %rightLeftShift = a1/(1-a1);RightShift = -a2/(1-a2);if (LeftShift>0 & RightShift>0)PeakShift = RightShift;elsePeakShift = LeftShift;endelseif(method==9) %Quinn method for the Hanning windowa1 = real(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %lefta2 = real(SpecV ec(PeakIdx+1)/SpecVec(PeakIdx)); %rightLeftShift = (2*a1+1)/(1-a1);RightShift = -(2*a2+1)/(1-a2);if (LeftShift>0 & RightShift>0)PeakShift = RightShift;elsePeakShift = LeftShift;endendend。