多种频谱校正方法及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 M
V 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 ωωωωωω++-???=+????
-?+???=+???+???
2121
cos()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 M
X X M m X X δ+++=+--6.加矩形窗情形,采用解析单频模型的复合复比值校正[3]
11Re ()M m M M X M m X X δ++??=+- ?-??
11m R m m X X X δ++=-,111
1m m L m m m m X X X X X X δ---=-=--0.5)0.5)m L m R
δδδδδ=-++（（7.加汉宁窗情形,采用解析单频模型的复合复比值校正[3]
1
12Re ()M M m M M X X M m X X δ++??+=+- ?-??
112m m R m m X X X X δ+++=-,1111
221m 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 R L R L R
ααδδαα++==---，, 00 R
R L R δδδδδ>>?=??当且其它
References
1.Schoukens,J.,R.Pintelon,H.Van Hamme.The interpolated fast Fourier transform:A
comparative 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.IEEE
Transactions on Signal Processing.1994,42(5):1264-1268.
%========================这是调用调试==================
DT=1;
N=1024;
PHI=pi/3;
Ampl=1;
CiR=11.9;%cycles in record
Freq=CiR/(DT*N);%frequency
TV=[0:N-1];
DatVec=Ampl*cos(Freq*TV*2*pi+PHI);
FV=fft(DatVec);
figure
subplot(2,1,1);plot(TV,DatVec);
subplot(2,1,2);plot(abs(FV(1:round(N/2.56))));grid on
[MV,MI]=max(abs(FV));
%加矩形窗的解析校正--1
FreqShift=SpecCorr(FV,MI,N,1);
%加矩形窗的解析单频模型校正--2
FreqShift=SpecCorr(FV,MI,N,2);
%加汉宁窗的解析单频模型校正--3
HanDat=DatVec.*hanning(N,'periodic')';
FV=fft(HanDat);
FreqShift=SpecCorr(FV,MI,N,3);
%加矩形窗的解析单频模型校正+复比值法--4
FV=fft(DatVec);
FreqShift=SpecCorr(FV,MI,N,4);
%加汉宁窗的解析单频模型校正+复比值法--5
HanDat=DatVec.*hanning(N,'periodic')';
FV=fft(HanDat);
FreqShift=SpecCorr(FV,MI,N,5);
%加矩形窗的解析单频模型校正+复比值法+左右平均--6
FV=fft(DatVec);
FreqShift=SpecCorr(FV,MI,N,6);
%加汉宁窗的解析单频模型校正+复比值法+左右平均--7
HanDat=DatVec.*hanning(N,'periodic')';
FV=fft(HanDat);
FreqShift=SpecCorr(FV,MI,N,7);
%加矩形窗的解析单频模型校正+Quinn算法--8
FV=fft(DatVec);
FreqShift=SpecCorr(FV,MI,N,8);
%加汉宁窗的解析单频模型校正+Quinn算法--9
HanDat=DatVec.*hanning(N,'periodic')';
FV=fft(HanDat);
FreqShift=SpecCorr(FV,MI,N,9);
===========这是子程序======================
%spectrum correction assemble
%the sampling interval is1s(or unitary)
%Input:SpecVec--Discrte Fourier Spectrum from FFT
%PeakIdx--the peak index,noting the matrix in MatLab start from1
%TL--the length(or the point number)of the FFT
%method--correction method
%output:PeakShift--the corrected peak shifting from the peak in discrete
%spectrum
function PeakShift=SpecCorr(SpecVec,PeakIdx,TL,method)
%picking up the second highest spectrum line
if(abs(SpecVec(PeakIdx-1))>abs(SpecVec(PeakIdx+1)))
IP=[PeakIdx-1,PeakIdx];
ShiftCorr=-1;%shift aligning with the PeakIdx
else
IP=[PeakIdx,PeakIdx+1];
ShiftCorr=0;%shift aligning with the PeakIdx
end
II=IP(1)-1;%noting that the index of a matrix in MATLAB starts from1,not zero if(method==1)%an analyitic solution for rectangular window
U=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 PeakIdx
elseif(method==3)%based on the analytical-single-tone model for Hanning window PeakShift=(2*abs(SpecVec(IP(2)))-abs(SpecVec(IP(1))))/(abs(SpecVec(IP(2)))+abs(SpecVec(IP(1) )));
PeakShift=PeakShift+ShiftCorr;%shift aligning with the PeakIdx
elseif(method==4)%based on the analytical-single-tone model for rectangular window with complex correction
PeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));
PeakShift=PeakShift+ShiftCorr;%shift aligning with the PeakIdx
elseif(method==5)%based on the analytical-single-tone model for Hanning window with complex correction
PeakShift=(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+average
PeakShift=real(SpecVec(IP(2))/(SpecVec(IP(2))-SpecVec(IP(1))));
MaxPeakShift=PeakShift+ShiftCorr;%shift aligning with the PeakIdx
LeftShift=real(SpecVec(PeakIdx)/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1)))-1;
RightShift=real(SpecVec(PeakIdx+1)/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx)));
%average
PeakShift=(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 PeakIdx LeftShift=(2*SpecVec(PeakIdx)+SpecVec(PeakIdx-1))/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1))-1;
RightShift=(2*SpecVec(PeakIdx+1)+SpecVec(PeakIdx))/(SpecVec(PeakIdx+1)-SpecVec(PeakIdx) );
%average
PeakShift=(0.5-MaxPeakShift)*LeftShift+(0.5+MaxPeakShift)*RightShift;
elseif(method==8)%Quinn method for the rectangular window
a1=real(SpecVec(PeakIdx-1)/SpecVec(PeakIdx));%left
a2=real(SpecVec(PeakIdx+1)/SpecVec(PeakIdx));%right
LeftShift=a1/(1-a1);
RightShift=-a2/(1-a2);
if(LeftShift>0&RightShift>0)
PeakShift=RightShift;
else
PeakShift=LeftShift;
end
elseif(method==9)%Quinn method for the Hanning window
a1=real(SpecVec(PeakIdx-1)/SpecVec(PeakIdx));%left
a2=real(SpecVec(PeakIdx+1)/SpecVec(PeakIdx));%right
LeftShift=(2*a1+1)/(1-a1);
RightShift=-(2*a2+1)/(1-a2);
if(LeftShift>0&RightShift>0)
PeakShift=RightShift;
else
PeakShift=LeftShift;
end
end
end
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 M V 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 ωωωωωω++-???=+?????-?+???=+???+??? 2121 cos()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] 1 1Re ()M M M X M m X X δ++??=+- ?-?? 5. 加汉宁窗情形,采用解析单频模型的复比值校正[3] 112()M M M M X X M m X X δ+++=+-- 6. 加矩形窗情形,采用解析单频模型的复合复比值校正[3]
用MATLAB进行FFT频谱分析
用MATLAB 进行FFT 频谱分析 假设一信号： ()()292.7/2cos 1.0996.2/2sin 1.06.0+++=t t R ππ 画出其频谱图。 分析： 首先，连续周期信号截断对频谱的影响。 DFT 变换频谱泄漏的根本原因是信号的截断。即时域加窗，对应为频域卷积，因此，窗函数的主瓣宽度等就会影响到频谱。 实验表明，连续周期信号截断时持续时间与信号周期呈整数倍关系时，利用DFT 变换可以得到精确的模拟信号频谱。举一个简单的例子： ()ππ2.0100cos +=t Y 其周期为0.02。截断时不同的持续时间影响如图一.1：（对应程序shiyan1ex1.m ） 图 错误！文档中没有指定样式的文字。.1 140.0160.0180.02 截断时，时间间期为周期整数倍，频谱图 0.0250.03 0100200300400500600 7008009001000 20 40 60 80 100 截断时，时间间期不为周期整数倍，频谱图
其次，采样频率的确定。 根据Shannon 采样定理，采样带限信号采样频率为截止频率的两倍以上，给定信号的采样频率应>1/7.92，取16。 再次，DFT 算法包括时域采样和频域采样两步，频域采样长度M 和时域采样长度N 的关系要符合M ≧N 时，从频谱X(k)才可完全重建原信号。 实验中信号R 经采样后的离散信号不是周期信号，但是它又是一个无限长的信号，因此处理时时域窗函数尽量取得宽一些已接近实际信号。 实验结果如图一.2：其中，0点位置的冲激项为直流分量0.6造成（对应程序为shiyan1.m ） 图 错误！文档中没有指定样式的文字。.2 ?ARMA （Auto Recursive Moving Average ）模型： 将平稳随机信号x(n)看作是零均值，方差为σu 2的白噪声u(n)经过线性非移变系统H(z)后的输出，模型的传递函数为 020406080100120140160180200 0.4 0.50.60.7 0.800.050.10.150.20.250.30.350.40.450.5 50100 150
计算方法_全主元消去法_matlab程序
%求四阶线性方程组的MA TLAB程序 clear Ab=[0.001 2 1 5 1; 3 - 4 0.1 -2 2; 2 -1 2 0.01 3; 1.1 6 2.3 9 4];%增广矩阵 num=[1 2 3 4];%未知量x的对应序号 for i=1:3 A=abs(Ab(i:4,i:4));%系数矩阵取绝对值 [r,c]=find(A==max(A(:))); r=r+i-1;%最大值对应行号 c=c+i-1;%最大值对应列号 q=Ab(r,:),Ab(r,:)=Ab(i,:),Ab(i,:)=q;%行变换 w=Ab(:,c),Ab(:,c)=Ab(:,i),Ab(:,i)=w;%列变换 n=num(i),num(i)=num(c),num(c)=n;%列变换引起未知量x次序变化for j=i:3 Ab(j+1,:)=-Ab(j+1,i)*Ab(i,:)/Ab(i,i)+Ab(j+1,:);%消去过程 end end %最后得到系数矩阵为上三角矩阵 %回代算法求解上三角形方程组 x(4)=Ab(4,5)/Ab(4,4); x(3)=(Ab(3,5)-Ab(3,4)*x(4))/Ab(3,3); x(2)=(Ab(2,5)-Ab(2,3)*x(3)-Ab(2,4)*x(4))/Ab(2,2); x(1)=(Ab(1,5)-Ab(1,2)*x(2)-Ab(1,3)*x(3)-Ab(1,4)*x(4))/Ab(1,1); for s=1:4 fprintf('未知量x%g =%g\n',num(s),x(s)) end %验证如下 %A=[0.001 2 1 5 1; 3 -4 0.1 -2 2;2 -1 2 0.01 3; 1.1 6 2.3 9 4]; %b=[1 2 3 4]'; %x=A\b; %x1= 1.0308 %x2= 0.3144 %x3= 0.6267 %x4= -0.0513
信号的频域分析及MATLAB实现.doc
《M A T L A B电子信息应用》 课程设计 设计五 信号的频域分析及MATLAB实现 学院： 专业： 班级： 姓名： 学号：
信号的频域分析及MATLAB实现 一、设计目的 通过该设计，理解傅里叶变换的定义及含义，掌握对信号进行频域分析的方法。 二、课程设计环境 计算机 MATLAB软件 三、设计内容及主要使用函数 快速傅里叶变换的应用 1）滤波器频率响应 对特定频率的频点或该频点以外的频率进行有效滤除的电路，就是滤波器。其功能就是得到一个特定频率或消除一个特定频率，滤波器是一种对信号有处理作用的器件或电路。主要作用是：让有用信号尽可能无衰减的通过，对无用信号尽可能大的。 滤波器的类型：巴特沃斯响应（最平坦响应），贝赛尔响应，切贝雪夫响应。 滤波器冲激响应的傅里叶变换就是该滤波器的频率响应。
2)快速卷积 卷积定理指出，函数卷积的傅里叶变换是函数傅里叶变换的乘积。即一个域中的卷积相当于另一个域中的乘积，例如时域中的卷积就对应于频域中的乘积。其中表示f 的傅里叶变换。 这一定理对拉普拉斯变换、双边拉普拉斯变换等各种傅里叶变换的变体同样成立。在调和分析中还可以推广到在局部紧致的阿贝尔群上定义的傅里叶变换。 利用卷积定理可以简化卷积的运算量。对于长度为n 的序列，按照卷积的定义进行计算，需要做2n - 1组对位乘法，其计算复杂度为；而利用傅里叶变换将序列变换到频域上后，只需要一组对位乘法，利用傅里叶变换的快速算法之后，总的计算复杂度为。这一结果可以在快速乘法计算中得到应用。 1. 信号的离散傅里叶变换 有限长序列的离散傅里叶变换公式为： kn N j N n e n x k X )/2(10)()(π--=∑= ∑==1_0)/2()(1)(N n kn N j e k X N n x π MATLAB 函数：fft 功能是实现快速傅里叶变换，fft 函数的格式为： ),(x fft y =返回向量x 的不连续fourier 变换。 若)6 cos()(πn n x =是一个N=12的有限序列，利用MATLAB 计算
离散频谱校正技术
图3.1.1 窗函数的频谱函数 三、离散频谱校正技术 经FFT 得到的离散频谱其幅值、相位和频率都可能产生较大的误差。从理论上分析，加矩形窗时单谐波频率的最大误差可达36.4％，即使加其它窗时，也不能完全消除此影响，如加Hanning 窗时，只进行幅值恢复时的最大误差仍高达15.3%，相位误差更大，高达90度。 目前国内外有四种对幅值谱或功率谱进行校正的方法：第一种方法是离散频谱能量重心校正法，第二种方法是对幅值谱进行校正的比值法，第三种方法是FFT+DFT 谱连续细化分析傅立叶变换法，第四种方法是相位差法，这些方法各有其特点。在相位差校正法中，有时移法、缩短窗长法和综合法。 1．比值校正法 这种方法利用频率归一化后差值为1的主瓣峰顶附近二条谱线的窗谱函数比值，建立一个以校正频率为变量的方程，解出校正频率，进而进行幅值和相位校正。解方程求校正频率的方法是多样化的，直接导出公式的方法称比值公式法，利用迭代求解的方法称为比值迭代公式法，用搜索求解的方法称比值峰值搜索法。研究表明，加Hanning 窗的比例校正法精度非常高，频率误差小于0.0001f ?，幅值误差小于万分之一，相位误差小于1度。 （1）频率校正 频率校正即求出主瓣中心的横坐标。设窗函数的频谱函数为 ()x f ，()x f 对称于y 轴，见图3.1.1。对于任一x ，窗谱函数为()x f ， 离散频谱为y x ；对于任一()1+x ，窗谱函数为()1+x f ，离散频谱为 y x +1，构造v 为间隔为1的两点()x f 、()1+x f 的比值函数，由()x f 、()1+x f 、y x 和y x +1就能求出x 。由于f(x)的函数表达式为已知，故 可构造一函数 v F x f x f x y y x x == +=+()() ()11 (3.1.1) v 是间隔为1的两点的比值，是x 的函数，对上式解出其反函数： x g v =() (3.1.2) 即求解谱线校正量x k x -=?=?，这种方法称为比值公式法。 校正频率为： N f k k f s x ) (?+= (3.1.3) 式中，()12/,,2,1,0-=N k k Λ为谱线号，N 为分析点数，s f 为采样频率。 （2）幅值校正 设窗函数的频谱模函数为()x f ，主瓣函数为： )(0x x Af y -= (3.1.4) 这就是信号频谱与窗函数卷积的结果，式中，A 为真实幅值，对应主瓣中心0x ，现将k y y =，k x =代入式(3.1.4)得： )(0x k Af y k -= (3.1.5)
用Matlab进行信号与系统的时、频域分析
课程实验报告 题目：用Matlab进行 信号与系统的时、频域分析 学院 学生姓名 班级学号 指导教师 开课学院 日期 用Matlab进行信号与系统的时、频域分析 一、实验目的 进一步了解并掌握Matlab软件的程序编写及运行； 掌握一些信号与系统的时、频域分析实例； 了解不同的实例分析方法，如：数值计算法、符号计算法； 通过使用不同的分析方法编写相应的Matlab程序； 通过上机，加深对信号与系统中的基本概念、基本理论和基本分析方法的理解。 二、实验任务 了解数值计算法编写程序，解决实例； 在Matlab上输入三道例题的程序代码，观察波形图； 通过上机实验，完成思考题； 完成实验报告。 三、主要仪器设备
硬件：微型计算机 软件：Matlab 四、 实验内容 （1） 连续时间信号的卷积 已知两个信号)2()1()(1---=t t t x εε和)1()()(2--=t t t x εε，试分别画出)(),(21t x t x 和卷积)()()(21t x t x t y *=的波形。 程序代码： T=0.01; t1=1;t2=2; t3=0;t4=1; t=0:T:t2+t4; x1=ones(size(t)).*((t>t1)-(t>t2)); x2=ones(size(t)).*((t>t3)-(t>t4)); y=conv(x1,x2)*T; subplot(3,1,1),plot(t,x1); ylabel('x1(t)'); subplot(3,1,2),plot(t,x2); ylabel('x2(t)'); subplot(3,1,3),plot(t,y(1:(t2+t4)/T+1)); ylabel('y(t)=x1*x2'); xlabel('----t/s'); （2）已知两个信号)()(t e t x t ε-=和)()(2/t te t h t ε-=，试用数值计算法求卷积，并分别画出)(),(t h t x 和卷积)()()(t h t x t y *=的波形。 程序代码： t2=3;t4=11; T=0.01; t=0:T:t2+t4; x=exp(-t).*((t>0)-(t>t2)); h=t.*exp(-t/2).*((t>0)-(t>t4)); y=conv(x,h)*T; yt=4*exp(-t)+2*t.*exp(-1/2*t)-4*exp(-1/2*t); subplot(3,1,1),plot(t,x); ylabel('x(t)'); subplot(3,1,2),plot(t,h); ylabel('h(t)'); subplot(3,1,3),plot(t,y(1:(t2+t4)/T+1),t,yt,'--r'); legend('by numberical','Theoretical'); ylabel('y=x*h'); xlabel('----t/s'); （3）求周期矩形脉冲信号的频谱图，已知s T s A 5.0,1.0,1===τ
多种频谱校正方法及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 M V 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 ωωωωωω++-???=+???? -?+???=+???+??? 2121 cos()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 M X X M m X X δ+++=+--6.加矩形窗情形,采用解析单频模型的复合复比值校正[3]
Matlab对采样数据进行频谱分析
使用Matlab对采样数据进行频谱分析 1、采样数据导入Matlab 采样数据的导入至少有三种方法。 第一就是手动将数据整理成Matlab支持的格式，这种方法仅适用于数据量比较小的采样。 第二种方法是使用Matlab的可视化交互操作，具体操作步骤为：File --> Import Data，然后在弹出的对话框中找到保存采样数据的文件，根据提示一步一步即可将数据导入。这种方法适合于数据量较大，但又不是太大的数据。据本人经验，当数据大于15万对之后，读入速度就会显著变慢，出现假死而失败。 第三种方法，使用文件读入命令。数据文件读入命令有textread、fscanf、load 等，如果采样数据保存在txt文件中，则推荐使用 textread命令。如 [a,b]=textread('data.txt','%f%*f%f'); 这条命令将data.txt中保存的数据三个三个分组，将每组的第一个数据送给列向量a，第三个数送给列向量b，第二个数据丢弃。命令类似于C语言，详细可查看其帮助文件。文件读入命令录入采样数据可以处理任意大小的数据量，且录入速度相当快，一百多万的数据不到20秒即可录入。强烈推荐！ 2、对采样数据进行频谱分析 频谱分析自然要使用快速傅里叶变换FFT了，对应的命令即 fft ，简单使用方法为：Y=fft(b,N)，其中b即是采样数据，N为fft数据采样个数。一般不指定N，即简化为Y=fft(b)。Y即为FFT变换后得到的结果，与b的元素数相等，为复数。以频率为横坐标，Y数组每个元素的幅值为纵坐标，画图即得数据b的幅频特性；以频率为横坐标，Y数组每个元素的角度为纵坐标，画图即得数据b的相频特性。典型频谱分析M程序举例如下： clc fs=100; t=[0:1/fs:100]; N=length(t)-1;%减1使N为偶数 %频率分辨率F=1/t=fs/N p=1.3*sin(0.48*2*pi*t)+2.1*sin(0.52*2*pi*t)+1.1*sin(0.53*2*pi*t)... +0.5*sin(1.8*2*pi*t)+0.9*sin(2.2*2*pi*t); %上面模拟对信号进行采样，得到采样数据p，下面对p进行频谱分析 figure(1) plot(t,p); grid on title('信号 p(t)'); xlabel('t') ylabel('p')
讲座1-3 离散频谱校正技术(DOC)
图3.1.1 窗函数的频谱函数 讲座1-3 三、离散频谱校正技术 经FFT 得到的离散频谱其幅值、相位和频率都可能产生较大的误差。从理论上分析，加矩形窗时单谐波频率的最大误差可达36.4％，即使加其它窗时，也不能完全消除此影响，如加Hanning 窗时，只进行幅值恢复时的最大误差仍高达15.3%，相位误差更大，高达90度。 目前国内外有四种对幅值谱或功率谱进行校正的方法：第一种方法是离散频谱能量重心校正法，第二种方法是对幅值谱进行校正的比值法，第三种方法是FFT+DFT 谱连续细化分析傅立叶变换法，第四种方法是相位差法，这些方法各有其特点。在相位差校正法中，有时移法、缩短窗长法和综合法。 1．比值校正法 这种方法利用频率归一化后差值为1的主瓣峰顶附近二条谱线的窗谱函数比值，建立一个以校正频率为变量的方程，解出校正频率，进而进行幅值和相位校正。解方程求校正频率的方法是多样化的，直接导出公式的方法称比值公式法，利用迭代求解的方法称为比值迭代公式法，用搜索求解的方法称比值峰值搜索法。研究表明，加Hanning 窗的比例校正法精度非常高，频率误差小于0.0001f ?，幅值误差小于万分之一，相位误差小于1度。 （1）频率校正 频率校正即求出主瓣中心的横坐标。设窗函数的频谱函数为 ()x f ，()x f 对称于y 轴，见图3.1.1。对于任一x ，窗谱函数为()x f ， 离散频谱为y x ；对于任一()1+x ，窗谱函数为()1+x f ，离散频谱为y x +1，构造v 为间隔为1的两点()x f 、()1+x f 的比值函数，由()x f 、 ()1+x f 、y x 和y x +1就能求出x 。由于f(x)的函数表达式为已知，故 可构造一函数 v F x f x f x y y x x == +=+()() ()11 (3.1.1) v 是间隔为1的两点的比值，是x 的函数，对上式解出其反函数： x g v =() (3.1.2) 即求解谱线校正量x k x -=?=?，这种方法称为比值公式法。 校正频率为： N f k k f s x ) (?+= (3.1.3) 式中，()12/,,2,1,0-=N k k 为谱线号，N 为分析点数，s f 为采样频率。 （2）幅值校正 设窗函数的频谱模函数为()x f ，主瓣函数为： )(0x x Af y -= (3.1.4) 这就是信号频谱与窗函数卷积的结果，式中，A 为真实幅值，对应主瓣中心0x ，现将k y y =，k x =代 入式(3.1.4)得：
(整理)matlab16常用计算方法.
常用计算方法 1．超越方程的求解 一超越方程为 x (2ln x – 3) -100 = 0 求超越方程的解。 [算法]方法一：用迭代算法。将方程改为 01002ln()3 x x =- 其中x 0是一个初始值，由此计算终值x 。取最大误差为e = 10-4，当| x - x 0| > e 时，就用x 的值换成x 0的值，重新进行计算；否则| x - x 0| < e 为止。 [程序]P1_1abs.m 如下。 %超越方程的迭代算法 clear %清除变量 x0=30; %初始值 xx=[]; %空向量 while 1 %无限循环 x=100/(2*log(x0)-3); %迭代运算 xx=[xx,x]; %连接结果 if length(xx)>1000,break ,end %如果项数太多则退出循环（暗示发散） if abs(x0-x)<1e-4,break ,end %当精度足够高时退出循环 x0=x; %替换初值 end %结束循环 figure %创建图形窗口 plot(xx,'.-','LineWidth',2,'MarkerSize',12)%画迭代线'.-'表示每个点用.来表示，再用线连接 grid on %加网格 fs=16; %字体大小 title('超越方程的迭代折线','fontsize',fs)%标题 xlabel('\itn','fontsize',fs) %x 标签 ylabel('\itx','fontsize',fs) %y 标签 text(length(xx),xx(end),num2str(xx(end)),'fontsize',fs)%显示结果 [图示]用下标作为自变量画迭代的折线。如P0_20_1图所示，当最大误差为10-4时，需要迭代19次才能达到精度，超越方程的解为27.539。 [算法]方法二：用求零函数和求解函数。将方程改为函数 100()2ln()3f x x x =-- MATLAB 求零函数为fzero ，fzero 函数的格式之一是 x = fzero(f,x0) 其中，f 表示求解的函数文件，x0是估计值。fzero 函数的格式之二是 x = fzero(f,[x1,x2])
matlab频谱分析
设计出一套完整的系统，对信号进行频谱分析和滤波处理； 1．产生一个连续信号，包含低频，中频，高频分量，对其进行采样，进行频谱分析，分别设计三种高通，低通，带通滤波器对信号进行滤波处理，观察滤波后信号的频谱。 2.采集一段含有噪音的语音信号（可以录制含有噪音的信号，或者录制语音后再加进噪音信号），对其进行采样和频谱分析，根据分析结果设计出一合适的滤波器滤除噪音信号。 %写上标题 %设计低通滤波器： [N,Wc]=buttord() %估算得到Butterworth低通滤波器的最小阶数N和3dB截止频率Wc [a,b]=butter(N,Wc); %设计Butterworth低通滤波器 [h,f]=freqz(); %求数字低通滤波器的频率响应 figure(2); % 打开窗口2 subplot(221); %图形显示分割窗口 plot(f,abs(h)); %绘制Butterworth低通滤波器的幅频响应图 title(巴氏低通滤波器''); grid; %绘制带网格的图像 sf=filter(a,b,s); %叠加函数S经过低通滤波器以后的新函数 subplot(222); plot(t,sf); %绘制叠加函数S经过低通滤波器以后的时域图形 xlabel('时间(seconds)'); ylabel('时间按幅度'); SF=fft(sf,256); %对叠加函数S经过低通滤波器以后的新函数进行256点的基—2快速傅立叶变换 w= %新信号角频率 subplot(223); plot()); %绘制叠加函数S经过低通滤波器以后的频谱图 title('低通滤波后的频谱图'); %设计高通滤波器 [N,Wc]=buttord() %估算得到Butterworth高通滤波器的最小阶数N和3dB截止频率Wc [a,b]=butter(N,Wc,'high'); %设计Butterworth高通滤波器 [h,f]=freqz(); %求数字高通滤波器的频率响应 figure(3); subplot(221); plot()); %绘制Butterworth高通滤波器的幅频响应图 title('巴氏高通滤波器'); grid; %绘制带网格的图像 sf=filter(); %叠加函数S经过高通滤波器以后的新函数 subplot(222); plot(t,sf); ;%绘制叠加函数S经过高通滤波器以后的时域图形 xlabel('Time(seconds)'); ylabel('Time waveform'); w; %新信号角频率 subplot(223);
matlab用于计算方法的源程序
1、Newdon迭代法求解非线性方程 function [x k t]=NewdonToEquation(f,df,x0,eps) %牛顿迭代法解线性方程 %[x k t]=NewdonToEquation(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:原函数，定义为内联函数 ?:函数的倒数，定义为内联函数 %x0:初始值 %eps:误差限 % %应用举例： %f=inline('x^3+4*x^2-10'); ?=inline('3*x^2+8*x'); %x=NewdonToEquation(f,df,1,0.5e-6) %[x k]=NewdonToEquation(f,df,1,0.5e-6) %[x k t]=NewdonToEquation(f,df,1,0.5e-6) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquation(f,df,1) if nargin==3 eps="0".5e-6; end tic; k=0; while 1 x="x0-f"(x0)./df(x0); k="k"+1; if abs(x-x0) < eps || k >30 break; end x0=x; end t=toc; if k >= 30 disp('迭代次数太多。'); x="0"; t="0"; end
2、Newdon迭代法求解非线性方程组 function y="NewdonF"(x) %牛顿迭代法解非线性方程组的测试函数 %定义是必须定义为列向量 y(1,1)=x(1).^2-10*x(1)+x(2).^2+8; y(2,1)=x(1).*x(2).^2+x(1)-10*x(2)+8; return; function y="NewdonDF"(x) %牛顿迭代法解非线性方程组的测试函数的导数 y(1,1)=2*x(1)-10; y(1,2)=2*x(2); y(2,1)=x(2).^+1; y(2,2)=2*x(1).*x(2)-10; return; 以上两个函数仅供下面程序的测试 function [x k t]=NewdonToEquations(f,df,x0,eps) %牛顿迭代法解非线性方程组 %[x k t]=NewdonToEquations(f,df,x0,eps) %x:近似解 %k:迭代次数 %t:运算时间 %f:方程组（事先定义） ?:方程组的导数（事先定义） %x0:初始值 %eps:误差限 % %说明：由于虚参f和df的类型都是函数，使用前需要事先在当前目录下采用函数M文件定义% 另外在使用此函数求解非线性方程组时，需要在函数名前加符号“@”，如下所示 % %应用举例： %x0=[0,0];eps=0.5e-6; %x=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps) %函数的最后一个参数也可以不写。默认情况下，eps=0.5e-6 %[x k t]=NewdonToEquations(@NewdonF,@NewdonDF,x0,eps)
MATLAB关于FFT频谱分析的程序
MATLAB关于FFT频谱分析的程序 %***************1.正弦波****************% fs=100;%设定采样频率 N=128; n=0:N-1; t=n/fs; f0=10;%设定正弦信号频率 %生成正弦信号 x=sin(2*pi*f0*t); figure(1); subplot(231); plot(t,x);%作正弦信号的时域波形 xlabel('t'); ylabel('y'); title('正弦信号y=2*pi*10t时域波形'); grid; %进行FFT变换并做频谱图 y=fft(x,N);%进行fft变换 mag=abs(y);%求幅值 f=(0:length(y)-1)'*fs/length(y);%进行对应的频率转换 figure(1); subplot(232); plot(f,mag);%做频谱图 axis([0,100,0,80]); xlabel('频率(Hz)'); ylabel('幅值');
title('正弦信号y=2*pi*10t幅频谱图N=128'); grid; %求均方根谱 sq=abs(y); figure(1); subplot(233); plot(f,sq); xlabel('频率(Hz)'); ylabel('均方根谱'); title('正弦信号y=2*pi*10t均方根谱'); grid; %求功率谱 power=sq.^2; figure(1); subplot(234); plot(f,power); xlabel('频率(Hz)'); ylabel('功率谱'); title('正弦信号y=2*pi*10t功率谱'); grid; %求对数谱 ln=log(sq); figure(1); subplot(235); plot(f,ln);
基于DSP的一种频谱校正算法的实现
第37卷第3期2014年6月电子器件Chinese Journal of Electron Devices Vol.37 No.3Jun.2014 收稿日期:2013-06-19 修改日期:2013-07-11A Spectrum Correction Algorithm Based on DSP TAN Chunhua ,YU Jian *,LIU Yahong ,QIAO Jinlong (Colleage of Computer and Control Engineering ,North University of China ,Taiyuan 030051,China )Abstract :Spectrum correction signal and information processing is important.The idea is to use an algorithm to search out a more precise spectrum peak,and its correction.This paper introduces the FFT /apFFT phase spectrum correction algorithm,built hardware platform based on TMS320C5535expounded spectrum correction algorithm for hardware and software implementation.This set uses a TLV320AIC3204A /D and D /A in one of the Codec chip and TMS320C5535internal DMA module signals occur in acquisition.Software using the TMS320C5535internal HWAFFT module 1024?point FFT computation,the CCS environment using C language programming of the FFT /apFFT phase spectrum correction algorithm,the program using the DMA data transfer,improved operational efficiency.The results show that the FFT /apFFT phase correction algorithm frequency spectrum estimation accuracy is high,the actual application of a certain reference value.Key words :frequency estimation;apFFT /FFT;HWAFFT;TLV320AIC3204EEACC :6140C doi :10.3969/j.issn.1005-9490.2014.03.041基于DSP 的一种频谱校正算法的实现 谭春花,禹 健*,刘亚翃,乔晋龙(中北大学计算机与控制工程学院,太原030051) 摘 要:频谱校正是信号与信息处理的重要内容三其思想是利用一种算法更精确的搜索出谱峰值,并对其进行校正三介绍了FFT /apFFT 相位差频谱校正算法原理,构建了基于TMS320C5535的硬件平台,阐述了频谱校正算法的硬件与软件实现三采用了TLV320AIC3204这一集A /D 与D /A 于一体的Codec 芯片以及TMS320C5535内部的DMA 模块实现信号的发生于采集三软件上利用TMS320C5535内部的HWAFFT 模块实现1024点的FFT 运算,在CCS 环境下利用C 语言编程实现了FFT /apFFT 相位差频谱校正算法,程序中利用DMA 进行数据的传送,提高了运行效率三运行结果表明了FFT /apFFT 相位差频谱校正算法频率估计精度高,对实际应用有一定的参考价值三 关键词:频谱校正;apFFT /FFT;HWAFFT;TLV320AIC3204 中图分类号:TP301.6 文献标识码:A 文章编号:1005-9490(2014)03-0570-04 在电力二铁路二通信二地质勘探二语音处理等应用 场合,存在着大量额关于信号的频率识别问题[1]三为估计频率的真实值,一种最直接的方法就是利用主谱线及其附近的几根旁谱线的幅值进行插值三Jain V K 首次提出基于矩形窗的插值方法[2],在此基础上,Grandke 应用了汉宁窗插值[3]三频谱校正方法的另外一个方向是基于相位差的比较法三文献[4]提出了利用分段FFT 的相位差提高正弦信号频率估计精度的方法三文献[5]中提出了一种时移相 位差校正法三文献[6]提出了基于DFT 相位的正弦 波频率的高精度估计方法三以上这些方法都是在传 统FFT 的架构下进行的,因此固有的频谱泄漏效应 会在很大程度上影响这些校正方法的精度三在文献[7]中王兆华教授提出了与传统FFT 谱分析相比,全相位FFT 谱分析具有更优良的抑制频谱泄漏的 性能三数字信号处理器DSP 是种可编程的高性能 处理器三文中利用TMS320C5535DSP 强大的数据 处理能力以及其特有的HWAFFT 模块,实现了 FFT /apFFT 相位差频谱校正算法,提高了运算效率,实际验证了该算法频率估计精度高三1 系统硬件结构系统设计以TI 公司的C55x 系列的一款C5535处理器为核心,辅助外围电路构成三DSP 负责对采集
Matlab频谱分析程序
Matlab 信号处理工具箱 谱估计专题 频谱分析 Spectral estimation （谱估计）的目标是基于一个有限的数据集合描述一个信号的功率（在频率上的）分布。功率谱估计在很多场合下都是有用的，包括对宽带噪声湮没下的信号的检测。 从数学上看，一个平稳随机过程n x 的power spectrum （功率谱）和correlation sequence （相关序列）通过discrete-time Fourier transform （离散时间傅立叶变换）构成联系。从normalized frequency （归一化角频率）角度看，有下式 ()()j m xx xx m S R m e ωω∞ -=-∞ = ∑ 注：()() 2 xx S X ωω=，其中( )/2 /2 lim N j n n N n N X x e ωω=-=∑ πωπ-<≤。其matlab 近似为X=fft(x,N)/sqrt(N)，在下文中()L X f 就是指matlab fft 函数的计算结果了 使用关系2/s f f ωπ=可以写成物理频率f 的函数，其中s f 是采样频率 ()()2/s jfm f xx xx m S f R m e π∞ -=-∞ = ∑ 相关序列可以从功率谱用IDFT 变换求得： ()()()/2 2//2 2s s s f jfm f j m xx xx xx s f S e S f e R m d df f πωπ π ωωπ--= =? ? 序列n x 在整个Nyquist 间隔上的平均功率可以表示为 ()()() /2 /2 02s s f xx xx xx s f S S f R d df f π π ωωπ--= =? ?
谱校正方法
谱校正方法 采样间隔归一化成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 M V 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 ωωωωωω++-???=+?????-?+???=+???+??? 2121 cos()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] 1 1Re ()M M M X M m X X δ++??=+- ?-?? 5. 加汉宁窗情形,采用解析单频模型的复比值校正[3] 112()M M M M X X M m X X δ+++=+-- 6. 加矩形窗情形,采用解析单频模型的复合复比值校正[3]
计算方法上机实验报告-MATLAB
《计算方法》实验报告 指导教师： 学院： 班级： 团队成员：
一、题目 例2.7应用Newton 迭代法求方程210x x --=在1x =附近的数值解 k x ，并使其满足8110k k x x ---< 原理： 在方程()0f x =解的隔离区间[],a b 上选取合适的迭代初值0x ，过曲线()y f x =的点()() 00x f x ，引切线 ()()()1000:'l y f x f x x x =+- 其与x 轴相交于点：()() 0100 'f x x x f x =-，进一步,过曲线()y f x =的 点()()11x f x , 引切线 ()()()2111: 'l y f x f x x x =+- 其与x 轴相交于点：() () 1211 'f x x x f x =- 如此循环往复，可得一列逼近方程()0f x =精确解*x 的点 01k x x x ,,,,,其一般表达式为： ()() 111 'k k k k f x x x f x ---=- 该公式所表述的求解方法称为Newton 迭代法或切线法。
程序： function y=f(x)%定义原函数 y=x^3-x-1; end function y1=f1(x0)%求导函数在x0点的值 syms x; t=diff(f(x),x); y1=subs(t,x,x0); end function newton_iteration(x0,tol)%输入初始迭代点x0及精度tol x1=x0-f(x0)/f1(x0);k=1;%调用f函数和f1函数 while abs(x1-x0)>=tol x0=x1;x1=x0-f(x0)/f1(x0);k=k+1; end fprintf('满足精度要求的数值为x(%d)=%1.16g\n',k,x1); fprintf('迭代次数为k=%d\n',k); end 结果：
利用matlab怎样进行频谱分析、、
利用matlab怎样进行频谱分析 图像的频率是表征图像中灰度变化剧烈程度的指标，是灰度在平面空间上的梯度。如：大面积的沙漠在图像中是一片灰度变化缓慢的区域，对应的频率值很低；而对于地表属性变换剧烈的边缘区域在图像中是一片灰度变化剧烈的区域，对应的频率值较高。傅立叶变换在实际中有非常明显的物理意义，设f是一个能量有限的模拟信号，则其傅立叶变换就表示f的谱。从纯粹的数学意义上看，傅立叶变换是将一个函数转换为一系列周期函数来处理的。从物理效果看，傅立叶变换是将图像从空间域转换到频率域，其逆变换是将图像从频率域转换到空间域。换句话说，傅立叶变换的物理意义是将图像的灰度分布函数变换为图像的频率分布函数，傅立叶逆变换是将图像的频率分布函数变换为灰度分布函数。 这样通过观察傅立叶变换后的频谱图，也叫功率图，我们首先就可以看出，图像的能量分布，如果频谱图中暗的点数更多，那么实际图像是比较柔和的（因为各点与邻域差异都不大，梯度相对较小），反之，如果频谱图中亮的点数多，那么实际图像一定是尖锐的，边界分明且边界两边像素差异较大的。对频谱移频到原点以后，可以看出图像的频率分布是以原点为圆心，对称分布的。将频谱移频到圆心除了可以清晰地看出图像频率分布以外，还有一个好处，它可以分离出有周期性规律的干扰信号，比如正弦干扰，一副带有正弦干扰，移频到原点的频谱图上可以看出除了中心以外还存在以某一点为中心，对称分布的亮点集合，这个集合就是干扰噪音产生的，这时可以很直观的通过在该位置放置带阻滤波器消除干扰。另外我还想说明以下几点： 1、图像经过二维傅立叶变换后，其变换系数矩阵表明： 若变换矩阵Fn原点设在中心，其频谱能量集中分布在变换系数短阵的中心附近(图中阴影区)。若所用的二维傅立叶变换矩阵Fn的原点设在左上角，那么图像信号能量将集中在系数矩阵的四个角上。这是由二维傅立叶变换本身性质决定的。同时也表明一股图像能量集中低频区域。 2、变换之后的图像在原点平移之前四角是低频，最亮，平移之后中间部分是低频，最亮，亮度大说明低频的能量大（幅角比较大）。 从计算机处理精度上就不难理解，一个长度为N的信号，最多只能有N/2+1个不同频率，再多的频率就超过了计算机所能所处理的精度范围）X[]数组又分两种，一种是表示余弦波的不同频率幅度值：Re X[]，另一种是表示正弦波的不同频率幅度值：Im X[]，Re是实数(Real)的意思，Im是虚数(Imagine)的意思，采用复数的表示方法把正余弦波组合起来进行表示，但这里我们不考虑复数的其它作用，只记住是一种组合方法而已，目的是为了便于表达（在后面我们会知道，复数形式的傅立叶变换长度是N，而不是N/2+1）。