9种谱校正方法及matlab代码

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 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]

11Re ()M m M M X M m X X δ++??=+- ?-?? 11m R m m

X 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 m

X 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 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,DatV ec);

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 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

% 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 from 1, 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(SpecV ec(IP(1))))/(abs(SpecV ec(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)+SpecV ec(PeakIdx-1))/(SpecVec(PeakIdx)-SpecVec(PeakIdx-1))-1;

RightShift=(2*SpecVec(PeakIdx+1)+SpecV ec(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(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %left

a2 = real(SpecV ec(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(SpecV ec(PeakIdx-1)/SpecVec(PeakIdx)); %left

a2 = real(SpecV ec(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作业二

6-1 试验目的计算特征值，实现算法 试验容：随机产生一个10阶整数矩阵，各数均在-5和5之间。 （1） 用MATLAB 函数“eig ”求矩阵全部特征值。 （2） 用幂法求A 的主特征值及对应的特征向量。 （3） 用基本QR 算法求全部特征值（可用MATLAB 函数“qr ”实现矩阵的QR 分解）。 原理 幂法：设矩阵A 的特征值为12n ||>||||λλλ≥???≥并设A 有完全的特征向量系12,,,n χχχ???(它们线性无关)，则对任意一个非零向量0n V R ∈所构造的向量序列1k k V AV -=有11()lim ()k j k k j V V λ→∞ -=， 其中()k j V 表示向量的第j 个分量。 为避免逐次迭代向量k V 不为零的分量变得很大（1||1λ>时）或很小（1||1λ<时），将每一步的k V 按其模最大的元素进行归一化。具体过程如下： 选择初始向量0V ，令1max(),,,1k k k k k k k V m V U V AU k m +===≥，当k 充分大时1111,max()max() k k U V χλχ+≈ ≈。 QR 法求全部特征值： 111 11222 111 ,1,2,3,k k k k k A A Q R R Q A Q R k R Q A Q R +++==????==??=???? ??????==?? 由于此题的矩阵是10阶的，上述算法计算时间过长，考虑采用改进算法——移位加速。迭 代格式如下： 1 k k k k k k k k A q I Q R A R Q q I +-=?? =+? 计算k A 右下角的二阶矩阵() () 1,1 1,() (),1 ,k k n n n n k k n n n n a a a a ----?? ? ??? 的特征值()()1,k k n n λλ-，当()()1,k k n n λλ-为实数时，选k q 为()()1,k k n n λλ-中最接近(),k n n a 的。 程序

计算方法_全主元消去法_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

频谱分析报告仪地使用方法

频谱分析仪的使用方法 13MHz信号。一般情况下，可以用示波器判断13MHz电路信号的存在与否，以及信号的幅度是否正常，然而，却无法利用示波器确定13MHz电路信号的频率是否正常，用频率计可以确定13MHz电路信号的有无，以及信号的频率是否准确，但却无法用频率计判断信号的幅度是否正常。然而，使用频谱分析仪可迎刃而解，因为频谱分析仪既可检查信号的有无，又可判断信号的频率是否准确，还可以判断信号的幅度是否正常。同时它还可以判断信号，特别是VCO信号是否纯净。可见频谱分析仪在手机维修过程中是十分重要的。 另外，数字手机的接收机、发射机电路在待机状态下是间隙工作的，所以在待机状态下，频率计很难测到射频电路中的信号，对于这一点，应用频谱分析仪不难做到。 一、使用前须知 在使用频谱分析仪之前，有必要了解一下分贝(dB)和分贝毫瓦(dBm)的基本概念，下面作一简要介绍。 1．分贝(dB) 分贝是增益的一种电量单位，常用来表示放大器的放大能力、衰减量等，表示的是一个相对量，分贝对功率、电压、电流的定义如下： 分贝数：101g(dB) 分贝数=201g(dB) 分贝数=201g(dB) 例如：A功率比B功率大一倍，那么，101gA／B=10182’3dB，也就是说，A功率比B功率大3dB， 2．分贝毫瓦(dBm) 分贝毫瓦(dBm)是一个表示功率绝对值的单位，计算公式为： 分贝毫瓦=101g(dBm) 例如，如果发射功率为lmw，则按dBm进行折算后应为：101glmw／1mw=0dBm。如果发射功率为40mw，则10g40w／1mw--46dBm。 二、频谱分析仪介绍 生产频谱分析仪的厂家不多。我们通常所知的频谱分析仪有惠普(现在惠普的测试设备分离出来，为安捷伦)、马可尼、惠美以及国产的安泰信。相比之下，惠普的频谱分析仪性能最好，但其价格也相当可观，早期惠美的5010频谱分析仪比较便宜，国产的安泰5010频谱分析仪的功能与惠美的5010差不多，其价格却便宜得多。 下面以国产安泰5010频谱分析仪为例进行介绍。 1．性能特点 AT5010最低能测到2.24uv，即是-100dBm。一般示波器在lmv，频率计要在20mv以上，跟频谱仪比相差10000倍。如用频率计测频率时，有的频率点测量很难，有的频率点测最不准，频率数字显示不稳定，甚至测不出来。这主要足频率计灵敏度问题，即信号低于20mv频率计就无能为力了，如用示波器测量时，信号5％失真示波器看不出来，在频谱仪上万分之一的失真都能看出来。

基于频域的校正方法及实验设计

2016届毕业（设计）论文 题目基于频域的校正方法及实验设计 专业班级过程自动化 学号 1204160134 学生姓名于春明 第一指导教师陈杰 指导教师职称 学院名称电气信息院 完成日期： 2016年 5月 20日

基于频域的校正方法及实验设计 Correction Method and Experimental Design based on Frequency Domain 学生姓名于春明 第一指导教师陈杰

摘要 在经典控制理论中，系统校正设计，就是在给定的性能指标下，对于给定的对象模型，确定一个能够完成系统满足的静态与动态性能指标要求的控制器，即确定校正器的结构与参数。串联校正控制器的频域设计方法中，使用的校正器有超前校正器、滞后校正器、滞后-超前校正器等。超前校正设计方法的特点是校正后系统的截止频率比校正前的大，系统的快速性能得到提高，这种校正设计方法对于要求稳定性好、超调量小以及动态过程响应快的系统被经常采用。滞后校正设计方法的特点是校正后系统的截止频率比校正前的小，系统的快速性能变差，但系统的稳定性能却得到提高，因此，在系统快速性要求不是很高，而稳定性与稳态精度要求很高的场合，滞后校正设计方法比较适合。滞后-超前校正设计是指既有滞后校正作用又有超前校正作用的校正器设计。它既具有了滞后校正高稳定性能、高精确度的好处，又具有超前校正响应快、超调小的优点，这种设计方法在要求较高的场合经常被采用。 关键词：

ABSTRACT In classical control theory, system design correction, that is, given the performance for a given object model to determine a controller to complete the system meets the static and dynamic performance requirements, namely to determine the structure and corrector parameter. Frequency Domain controller series corrected using correction has lead corrector corrector lag lag - lead Corrector like. Features advanced design correction method is to correct the system cut-off frequency than the fast performance is improved before the correction is large, the system, this correction method for design requires good stability, small overshoot and dynamic process of rapid response systems are often use. Lag compensation design approach is characterized by system after correcting the cutoff frequency than the fast performance before correction is small, the system is deteriorated but stable performance of the system has improved, therefore, the system speed requirements are not high, and the stability and steady-state high precision of the occasion, lag compensation design method is more suitable. Lag - lead Corrector design means there are both lagging corrective action ahead of corrective action calibration designs. It has a lag correction high stability, high accuracy advantages, but also has the leading correction fast response, small overshoot of the advantages of this design approach in demanding situations are often used. Keywords：

离散频谱校正技术

图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实现课后代码

第一章插值方法 1.1Lagrange插值 1.2逐步插值 1.3分段三次Hermite插值 1.4分段三次样条插值 第二章数值积分 2.1 Simpson公式 2.2 变步长梯形法 2.3 Romberg加速算法 2.4 三点Gauss公式 第三章常微分方程德差分方法 3.1 改进的Euler方法 3.2 四阶Runge-Kutta方法 3.3 二阶Adams预报校正系统 3.4 改进的四阶Adams预报校正系统 第四章方程求根 4.1 二分法 4.2 开方法 4.3 Newton下山法 4.4 快速弦截法 第五章线性方程组的迭代法 5.1 Jacobi迭代 5.2 Gauss-Seidel迭代 5.3 超松弛迭代 5.4 对称超松弛迭代 第六章线性方程组的直接法 6.1 追赶法 6.2 Cholesky方法 6.3 矩阵分解方法 6.4 Gauss列主元消去法

第一章插值方法 1.1Lagrange插值 计算Lagrange插值多项式在x=x0处的值. MATLAB文件：（文件名：Lagrange_eval.m）function [y0,N]= Lagrange_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Lagrange插值多项式在x0处的值 %N是Lagrange插值函数的权系数 m=length(X); N=zeros(m,1); y0=0; for i=1:m N(i)=1; for j=1:m if j~=i; N(i)=N(i)*(x0-X(j))/(X(i)-X(j)); end end y0=y0+Y(i)*N(i); end 用法》X=[…];Y=[…]; 》x0= ; 》[y0,N]= Lagrange_eval(X,Y,x0) 1.2逐步插值 计算逐步插值多项式在x=x0处的值. MATLAB文件：（文件名：Neville_eval.m）function y0=Neville_eval(X,Y,x0) %X,Y是已知插值点坐标 %x0是插值点 %y0是Neville逐步插值多项式在x0处的值 m=length(X); P=zeros(m,1); P1=zeros(m,1); P=Y; for i=1:m P1=P; k=1; for j=i+1:m k=k+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]

讲座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、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)

0计算方法及MATLAB实现简明讲义课件PPS8-1欧拉龙格法

第8章 常微分方程初值问题数值解法 8.1 引言 8.2 欧拉方法 8.3 龙格-库塔方法 8.4 单步法的收敛性与稳定性 8.5 线性多步法

8.1 引 言 考虑一阶常微分方程的初值问题 00(,),[,],(). y f x y x a b y x y '=∈=（1.1） （1.2） 如果存在实数 ，使得 121212(,)(,).,R f x y f x y L y y y y -≤-?∈（1.3） 则称 关于 满足李普希茨(Lipschitz ）条件， 称为 的李普希茨常数（简称Lips.常数）. 0>L f y L f (参阅教材386页)

计算方法及MATLAB 实现 所谓数值解法，就是寻求解 在一系列离散节点 )(x y <<<<<+121n n x x x x 上的近似值 . ,,,,,121+n n y y y y 相邻两个节点的间距 称为步长. n n n x x h -=+1 如不特别说明，总是假定 为定数， ),2,1( ==i h h i 这时节点为 . ) ,2,1,0(0 =+=i nh x x n 初值问题(1.1),(1.2)的数值解法的基本特点是采取 “步进式”. 即求解过程顺着节点排列的次序一步一步地向前推进. 00(,),[,], (). y f x y x a b y x y '=∈=

描述这类算法，只要给出用已知信息 ,,,21--n n n y y y 计算 的递推公式. 1+n y 一类是计算 时只用到前一点的值 ，称为单步法. 1+n y n y 另一类是用到 前面 点的值 ， 1+n y k 11,,,+--k n n n y y y 称为 步法. k 其次，要研究公式的局部截断误差和阶，数值解 与 精确解 的误差估计及收敛性，还有递推公式的计算 稳定性等问题. n y )(n x y 首先对方程 离散化，建立求数值解的递推 公式. ),(y x f y ='

谱校正方法

谱校正方法 采样间隔归一化成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实现结构可靠度计算.

用MATLAB实现结构可靠度计算 口徐华…朝泽刚‘u刘勇‘21 。 (【l】中国地质大学(武汉工程学院湖北?武汉430074; 12】河海大学土木工程学院江苏?南京210098 摘要:Matlab提供了各种矩阵的运算和操作,其中包含结构可靠度计算中常用的各种数值计算方法工具箱,本文从基本原理和相关算例分析两方面,阐述利用Matlab,编制了计算结构可靠度Matlab程.序,使得Matlab-语言在可靠度计算中得到应用。 关键词:结构可靠度Matlab软件最优化法 中图分类号:TP39文献标识码:A文章编号:1007-3973(200902-095-Ol 1结构可靠度的计算方法 当川概率描述结构的可靠性时,计算结构可靠度就是计算结构在规定时问内、规定条件F结构能够完成预定功能的概率。 从简单到复杂或精确稃度的不同,先后提出的可靠度计算方法有一次二阶矩方法、二次二阶矩方法、蒙特卡洛方法以及其他方法。一次■阶矩方法又分为。I-心点法和验算点法,其中验算点法足H前可靠度分析最常川的方法。 2最优化方法计算可靠度指标数学模型 由结构111n个任意分布的独立随机变量一,x:…以表示的结构极限状态方程为:Z=g(■.托…t=0,采用R-F将非正念变量当罱正态化,得到等效正态分布的均值o:和标准差虹及可靠度指标B,由可靠度指标B的几何意义知。o;辟

开始时验算点未知,把6看成极限状态曲面上点P(■,爿:---37,的函数,通过优化求解,找到B最小值。求解可靠皮指标aJ以归结为以下约束优化模型: rain睁喜t华,2 s.,.Z=g(工i,x2’,…,工:=0 如极限状态方栉巾某个变最(X。可用其他变量表示,则上述模型jfIJ‘转化为无约束优化模型: 。。B!:手f生丛r+阻:坚:坠:盐尘}二剐 t∞oY?’【叫,J 3用MATLAB实现结构可靠度计算 3.1Matlab简介 Matlab是++种功能强、效率高、便.丁.进行科学和工程计算的交互式软件包,汇集了人量数学、统计、科学和工程所需的函数,MATI.AB具有编程简甲直观、用户界mf友善、开放性强等特点。将MATLAB用于蒙特卡罗法的一个显著优点是它拥有功能强大的随机数发生器指令。 3.2算例 3.2.I例:已知非线形极限状态方程z=g(t r'H=567f r-0.5H2=0’f、r服从正态分布。IIf=0.6,o r=0.0786;la|_ 2.18,o r_0.0654;H服从对数正态分布。u H= 3218,O。 =0.984。f、r、H相互独立,求可靠度指标B及验算点(,,r’,H‘。 解:先将H当量正念化:h=ln H服从正态分布,且 ,‘-““了:等专虿’=,。49?口二-、『五ir面_。。3

频谱分析

频谱分析 利用傅里叶变换的方法对振动的信号进行分解，并按频率顺序展开，使其成为频率的函数，进而在频率域中对信号进行研究和处理的一种过程，称为频谱分析。 怎样进行频谱分析： 利用频谱分析仪进行测量，输入信号不能有失真，因此要按特定应用的要求设置频谱分析仪和优化测量步骤，以达到最好的技术指标。下面的测量提示对这些步骤有详细的说明。 1. 选择最好的分辨率带宽 (RBW) 必须认真考虑分辨率带宽 (RBW)的设置，因为他关系到频谱成分的分离，适宜的噪声基底的设置和信号的解调。 通过低电平信号的测量，可以看到使用窄RBW的优点。在使用窄RBW时，频谱分析仪显示出较低的平均噪声级 (DANL)，且动态范围增加，灵敏度有所改进。在图3中，把RBW从100kHz改变到10kHz 将能更好地分辨-95dBm的信号。 但并非任何情况都是最窄的RBW最好。对于调制信号，RBW一定要设置得足够宽，使它能将信号边带包括在内。如果忽略这一点，测量将是极不精确的。窄RBW设置的一项重要缺点是扫频速度。更宽的RBW设置在给定频率范围内允许更快的扫频。图4和图5比较了在200MHz频率范围内，10kHz和 3kHzRBW的扫频时间。 一定要知道RBW 选择时所必须的基本权衡因素，使得用户在明白哪些参数最为重要的时候，给以适当的优化。但在权衡不可避免时，

现代频谱分析仪可为您提供弱化，甚至消除这些因素的方法。通过使用数字信号处理，频谱分析仪在实现更精确的测量的同时还提供更高的速度，即使是使用窄RBW。 2. 改进测量精度 在进行任何测量前，必须了解有哪些可以改进幅度和频率测量精度的技术。 自校准功能可用来产生误差校正系数 (例如幅度改变—分辨率带宽)，分析仪随后用它校正测量数据，得到更好的幅度测量结果，并使您能在测量过程中更灵活地改变控制。 当被测装置接到经校准的分析仪时，信号传输网络可能会使感兴趣信号减弱或变形，必须在测量中排除这一影响，见图6。一种方法是使用分析仪的内置幅度校正功能，一个信号源以及一个功率表。图7给出了一个对DUT信号产生衰减的信号传输网络的频率响应。为消除这一有害效应，可在测量范围内若干存在问题的频率点上测量信号传输网络的衰减或增益。幅度校正给出频率—幅度表，用直线连接这些点得到“校正”波形，然后按这些校正值对输入信号进行偏置。在图8 中，信号传输网络不需要的衰减和增益已从测量中消除，

计算方法上机实验报告-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实现第二章作业

作者：夏云木子 1、 >> syms re(x) re(y) re(z) >> input('计算相对误差：'),re(x)=10/1991,re(y)=0.0001/1.991,re(y)=0.0000001/0.0001991 所以可知re(y)最小，即y精度最高 2、 >> format short,A=sqrt(2) >> format short e,B=sqrt(2) >> format short g,C=sqrt(2)

>> format long,D=sqrt(2) >> format long e,E=sqrt(2) >> format long g,F=sqrt(2) >> format bank,H=sqrt(2) >> format hex,I=sqrt(2) >> format +,J=sqrt(2) >> format,K=sqrt(2)

3、 >> syms A >> A=[sqrt(3) exp(7);sin(5) log(4)];vpa(pi*A,6) 4、1/6251-1/6252=1/6251*6252 5、（1）1/（1+3x）-(1-x)/(1+x)=x*(3*x-1)/[(1+3*x)*(1+x)] (2) sqrt(x+1/x)-sqrt(x-1/x)=2/x/[sqrt(x-1/x)+sqrt(x+1/x)] (3) log10(x1)-log(x2)=log10(x1/x2) (4) [1-cos(2*x)]/x =x^2/factorial(2)-x^4/factorial(4)+x^6/factorial(6)-…

数值计算方法matlab程序

function [x0,k]=bisect1(fun1,a,b,ep) if nargin<4 ep=1e-5; end fa=feval(fun1,a); fb=feval(fun1,b); if fa*fb>0 x0=[fa,fb]; k=0; return; end k=1; while abs(b-a)/2>ep x=(a+b)/2; fx=feval(fun1,x); if fx*fa<0 b=x; fb=fx; else a=x; fa=fx;

end end x0=(a+b)/2; >> fun1=inline('x^3-x-1'); >> [x0,k]=bisect1(fun1,1.3,1.4,1e-4) x0 = 1.3247 k = 7 >> 简单迭代法 function [x0,k]=iterate1(fun1,x0,ep,N) if nargin<4 N=500; end if nargin<3 ep=1e-5; end x=x0; x0=x+2*ep;

while abs(x-x0)>ep & k> fun1=inline('(x+1)^(1/3)'); >> [x0,k]=iterate1(fun1,1.5) x0 = 1.3247 k = 7 >> fun1=inline('x^3-1'); >> [x0,k]=iterate1(fun1,1.5) x0 = Inf k =