disp(' ')
disp(' srs.m ver 2.0 July 3, 2006')
disp(' by Tom Irvine Email: tomirvine@https://www.doczj.com/doc/7e13599335.html,')
disp(' ')
disp(' This program calculates the shock response spectrum')
disp(' of an acceleration time history, which is pre-loaded into Matlab.') disp(' The time history must have two columns: time(sec) & acceleration') disp(' ')
%
clear t;
clear y;
clear yy;
clear n;
clear fn;
clear a1;
clear a2
clear b1;
clear b2;
clear jnum;
clear THM;
clear resp;
clear x_pos;
clear x_neg;
%
iunit=input(' Enter acceleration unit: 1= G 2= m/sec^2 ');
%
disp(' ')
disp(' Select file input method ');
disp(' 1=external ASCII file ');
disp(' 2=file preloaded into Matlab ');
file_choice = input('');
%
if(file_choice==1)
[filename, pathname] = uigetfile('*.*');
filename = fullfile(pathname, filename);
%
fid = fopen(filename,'r');
THM = fscanf(fid,'%g %g',[2 inf]);
THM=THM';
else
THM = input(' Enter the matrix name: ');
end
%
t=double(THM(:,1));
y=double(THM(:,2));
%
tmx=max(t);
tmi=min(t);
n = length(y);
%
out1 = sprintf('\n %d samples \n',n);
disp(out1)
%
dt=(tmx-tmi)/(n-1);
sr=1./dt;
%
out1 = sprintf(' SR = %g samples/sec dt = %g sec \n',sr,dt); disp(out1)
%
fn(1)=input(' Enter the starting frequency (Hz) ');
if fn(1)>sr/30.
fn(1)=sr/30.;
end
%
idamp=input(' Enter damping format: 1= damping ratio 2= Q '); %
disp(' ')
if(idamp==1)
damp=input(' Enter damping ratio (typically 0.05) ');
else
Q=input(' Enter the amplification factor (typically Q=10) '); damp=1./(2.*Q);
end
%
disp(' ')
disp(' Select algorithm: ')
disp(' 1=Kelly-Richman 2=Smallwood ');
ialgorithm=input(' ');
%
tmax=(tmx-tmi) + 1./fn(1);
limit = round( tmax/dt );
n=limit;
yy=zeros(1,limit);
for i=1:length(y)
yy(i)=y(i);
end
%
disp(' ')
disp(' Calculating response..... ')
%
% SRS engine
%
for j=1:1000
%
omega=2.*pi*fn(j);
omegad=omega*sqrt(1.-(damp^2));
cosd=cos(omegad*dt);
sind=sin(omegad*dt);
domegadt=damp*omega*dt;
%
if(ialgorithm==1)
a1(j)=2.*exp(-domegadt)*cosd;
a2(j)=-exp(-2.*domegadt);
b1(j)=2.*domegadt;
b2(j)=omega*dt*exp(-domegadt);
b2(j)=b2(j)*( (omega/omegad)*(1.-2.*(damp^2))*sind -2.*damp*cosd );
b3(j)=0;
%
else
E=exp(-damp*omega*dt);
K=omegad*dt;
C=E*cos(K);
S=E*sin(K);
Sp=S/K;
%
a1(j)=2*C;
a2(j)=-E^2;
b1(j)=1.-Sp;
b2(j)=2.*(Sp-C);
b3(j)=E^2-Sp;
end
forward=[ b1(j), b2(j), b3(j) ];
back =[ 1, -a1(j), -a2(j) ];
%
resp=filter(forward,back,yy);
%
x_pos(j)= max(resp);
x_neg(j)= min(resp);
%
jnum=j;
if fn(j) > sr/8.
break
end
fn(j+1)=fn(1)*(2. ^ (j*(1./12.)));
end
%
% Output options
%
disp(' ')
disp(' Select output option ');
choice=input(' 1=plot only 2=plot & output text file ' );
disp(' ')
%
if choice == 2
%%
[writefname, writepname] = uiputfile('*','Save SRS data as');
writepfname = fullfile(writepname, writefname);
writedata = [fn' x_pos' (abs(x_neg))' ];
fid = fopen(writepfname,'w');
fprintf(fid,' %g %g %g\n',writedata');
fclose(fid);
%%
% disp(' Enter output filename ');
% SRS_filename = input(' ','s');
%
% fid = fopen(SRS_filename,'w');
% for j=1:jnum
% fprintf(fid,'%7.2f %10.3f %10.3f
\n',fn(j),x_pos(j),abs(x_neg(j)));
% end
% fclose(fid);
end
%
% Plot SRS
%
disp(' ')
disp(' Plotting output..... ')
%
% Find limits for plot
%
srs_max = max(x_pos);
if max( abs(x_neg) ) > srs_max
srs_max = max( abs(x_neg ));
end
srs_min = min(x_pos);
if min( abs(x_neg) ) < srs_min
srs_min = min( abs(x_neg ));
end
%
figure(1);
plot(fn,x_pos,fn,abs(x_neg),'-.');
%
if iunit==1
ylabel('Peak Accel (G)');
else
ylabel('Peak Accel (m/sec^2)');
end
xlabel('Natural Frequency (Hz)');
Q=1./(2.*damp);
out5 = sprintf(' Acceleration Shock Response Spectrum Q=%g ',Q);
title(out5);
grid;
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log ','YScale','log');
legend ('positive','negative',2);
%
ymax= 10^(round(log10(srs_max)+0.8));
ymin= 10^(round(log10(srs_min)-0.6));
%
fmax=max(fn);
fmin=fmax/10.;
%
fmax= 10^(round(log10(fmax)+0.5));
%
if fn(1) >= 0.1
fmin=0.1;
end
if fn(1) >= 1
fmin=1;
end
if fn(1) >= 10
fmin=10;
end
if fn(1) >= 100
fmin=100;
end
axis([fmin,fmax,ymin,ymax]);
%
disp(' ')
disp(' Plot pseudo velocity? ');
vchoice=input(' 1=yes 2=no ' );
if(vchoice==1)
figure(2);
%
% Convert to pseudo velocity
%
for j=1:jnum
if iunit==1
x_pos(j)=386.*x_pos(j)/(2.*pi*fn(j));
x_neg(j)=386.*x_neg(j)/(2.*pi*fn(j));
else
x_pos(j)=x_pos(j)/(2.*pi*fn(j));
x_neg(j)=x_neg(j)/(2.*pi*fn(j));
end
end
%
srs_max = max(x_pos);
if max( abs(x_neg) ) > srs_max
srs_max = max( abs(x_neg ));
end
srs_min = min(x_pos);
if min( abs(x_neg) ) < srs_min
srs_min = min( abs(x_neg ));
end
%
plot(fn,x_pos,fn,abs(x_neg),'-.');
%
if iunit==1
ylabel('Velocity (in/sec)');
else
ylabel('Velocity (m/sec)');
end
xlabel('Natural Frequency (Hz)');
Q=1./(2.*damp);
out5 = sprintf(' Pseudo Velocity Shock Response Spectrum Q=%g ',Q); title(out5);
grid;
set(gca,'MinorGridLineStyle','none','GridLineStyle',':','XScale','log ','YScale','log');
legend ('positive','negative',2);
%
ymax= 10^(round(log10(srs_max)+0.8));
ymin= 10^(round(log10(srs_min)-0.6));
%
axis([fmin,fmax,ymin,ymax]); end
功率谱估计及其MATLAB仿真 1经典功率谱估计 经典功率谱估计是将数据工作区外的未知数据假设为零,相当于数据加窗。经典功率谱估计方法分为:相关函数法(BT法)、周期图法以及两种改进的周期图估计法即平均周期图法和平滑平均周期图法,其中周期图法应用较多,具有代表性。 1.1相关函数法(BT法) 该方法先由序列x(n)估计出自相关函数R(n),然后对R(n)进行傅立叶变换,便得到x(n)的功率谱估计。当延迟与数据长度相比很小时,可以有良好的估计精度。 Matlab代码示例1(Btfangfa.M)： Fs=500;%采样频率 n=0:1/Fs:1; xn=cos(2*pi*40*n)+3*cos(2*pi*90*n)+randn(size(n));%产生含有噪声的序列 nfft=512; cxn=xcorr(xn,'unbiased');%计算序列的自相关函数 CXk=fft(cxn,nfft); Pxx=abs(CXk); index=0:round(nfft/2-1); %Round towards nearest integer. k=index*Fs/nfft; plot_Pxx=10*log10(Pxx(index+1)); figure(1); plot(k,plot_Pxx); 结果如下： 1.2周期图法(periodogram) 周期图法是把随机序列x(n)的N个观测数据视为一能量有限的序列,直接计算x(n)的离散傅立叶变换得X(k),然后再取其幅值的平方,并除以N,作为序列x(n)真实功率谱的估计。Matlab代码示例2(PEriod.M)： Fs=600; n=0:1/Fs:1;
冲击响应谱校准技术的研究 厉巍 陈永久 朱永晓 （贵州航天计量测试技术研究所，贵州 贵阳550009） 摘要：冲击响应谱试验已经成为大多数航天产品必做的力学环境试验项目之一，传统的冲击试验缺乏对冲击环境模拟的真实性，本文介绍了冲击响应谱的原理和冲击响应谱试验设备；用labVIEW 为平台，编写了冲击响应谱校准软件，为冲击响应谱试验机的校准与数据分析提供了通用性较好的校准分析方法，并基于PXI 系统设计了冲击响应谱校准装置。 关键词：航天产品LabVIEW 冲击响应谱 校准 PXI 系统 0引言 冲击响应谱试验机是用于完成冲击响应谱试验的环境试验设备，冲击响应谱是对产品实施抗冲击设计的分析基础，也是生产中冲击环境模拟试验的基本参数，在航空、航天重点型号科研生产及有关重大科技专项中，冲击响应谱试验已经成为必做的环境试验之一。产品在实际应用过程中受力情况复杂，其中，冲击激励会使设备激起强迫振动和固有频率响应，使产品性能和结构强度受到不同程度的损害甚至失效。航空、航天、电子等行业产品在生产、运输等过程中存在着各种冲击，而这对产品的质量和可靠性有着很大的负面影响。为了解决这一问题，在此基础上产生并发展起了冲击试验。近年来，随着对环境试验的认识不断提高，对冲击环境的模拟也提出了更高的要求，冲击响应谱试验也来越被关注。 1 冲击响应谱原理 冲击信号与一般的振动信号在许多方面具有不同的特性，工程中研究冲击信号的目的并不是研究冲击波形本身，而是更加注重冲击作用于系统的效果，或者说是研究冲击运动对系统的损伤势。不论用冲击的时间历程还是用频谱都难以描述冲击的损伤势，因此必须使用能够衡量冲击效果的冲击响应谱。 冲击响应谱系指一单自由度质量弹簧阻尼系统，当公共基础受到冲击激励时产生的响应峰值作为单自由度系统固有频率的函数绘出的图，其物理模型如图1所示。 图1 冲击响应谱的物理模型 数学模型可归结为如下微分方程的解： 式中，u x -=δ；
%求四阶线性方程组的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实现功率谱密度分析psd及详细解说 功率谱密度幅值的具体含义？？ 求信号功率谱时候用下面的不同方法，功率谱密度的幅值大小相差很大！ 我的问题是，计算具体信号时，到底应该以什么准则决定该选用什么方法啊？ 功率谱密度的幅植的具体意义是什么？？下面是一些不同方法计算同一信号的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);
I C S17.160 J04 中华人民共和国国家标准 G B/T29716.4 2018/I S O18431-4:2007 机械振动与冲击信号处理 第4部分:冲击响应谱分析 M e c h a n i c a l v i b r a t i o na n d s h o c k S i g n a l p r o c e s s i n g P a r t4:S h o c k-r e s p o n s e s p e c t r u ma n a l y s i s (I S O18431-4:2007,I D T) 2018-03-15发布2018-10-01实施中华人民共和国国家质量监督检验检疫总局
目 次 前言Ⅲ 引言Ⅳ 1 范围1 2 规范性引用文件1 3 术语和定义1 4 符号和缩略语1 5 冲击响应谱基本原理2 6 冲击响应谱的计算5 7 采样频率的影响9 参考文献12
前言 G B/T29716‘机械振动与冲击信号处理“由以下部分组成: 第1部分:引论; 第2部分:傅立叶变换的时域窗; 第3部分:时频分析方法; 第4部分:冲击响应谱分析; 第5部分:时基分析方法三 本部分为G B/T29716的第4部分三 本部分按照G B/T1.1 2009给出的规则起草三 本部分使用翻译法等同采用I S O18431-4:2007‘机械振动与冲击信号处理第4部分:冲击响应谱分析“三 与本部分中规范性引用的国际文件有一致性对应关系的我国文件如下: G B/T2298 2010机械振动二冲击与状态监测词汇(I S O2041:2009,I D T)三 本部分由全国机械振动二冲击与状态监测标准化技术委员会(S A C/T C53)提出并归口三 本部分起草单位:西北机电工程研究所二杭州亿恒科技有限公司二中国测试技术研究院二交通运输部公路科学研究所二孝感松林国际计测器有限公司二湖北省电力公司电力科学研究院二中船重工第七一一研究所三 本部分主要起草人:李超位二焦明纲二顾国富二王宝元二洪丽娜二赵玉刚三
功率谱估计性能分析及Matlab 仿真 1 引言 随机信号在时域上是无限长的，在测量样本上也是无穷多的，因此随机信号的能量是无限的，应该用功率信号来描述。然而，功率信号不满足傅里叶变换的狄里克雷绝对可积的条件，因此严格意义上随机信号的傅里叶变换是不存在的。因此，要实现随机信号的频域分析，不能简单从频谱的概念出发进行研究，而是功率谱[1]。 信号的功率谱密度描述随机信号的功率在频域随频率的分布。利用给定的 N 个样本数据估计一个平稳随机信号的功率谱密度叫做谱估计。谱估计方法分为两大类：经典谱估计和现代谱估计。经典功率谱估计如周期图法、自相关法等，其主要缺陷是描述功率谱波动的数字特征方差性能较差，频率分辨率低。方差性能差的原因是无法获得按功率谱密度定义中求均值和求极限的运算[2]。分辨率低的原因是在周期图法中，假定延迟窗以外的自相关函数全为0。这是不符合实际情况的，因而产生了较差的频率分辨率。而现代谱估计的目标都是旨在改善谱估计的分辨率，如自相关法和Burg 法等。 2 经典功率谱估计 经典功率谱估计是截取较长的数据链中的一段作为工作区，而工作区之外的数据假设为0，这样就相当将数据加一窗函数，根据截取的N 个样本数据估计出其功率谱[1]。 周期图法( Periodogram ) Schuster 首先提出周期图法。周期图法是根据各态历经的随机过程功率谱的定义进行的谱估计。 取平稳随机信号()x n 的有限个观察值(0),(1),...,(1)x x x n -，求出其傅里叶变换 1 ()()N j j n N n X e x n e ω ω---==∑ 然后进行谱估计
冲击响应谱 1简介 冲击响应谱通常简称“冲击谱”，它是工程中广泛应用的一个重要概念。国家电工委员会（IEC）、国家标准化组织(ISO)所属的技术委员会以及我国的国家标准，都已经把冲击谱作为规定冲击环境的方法之一。因此，冲击谱是对设备实施抗冲击设计的分析基础，也是控制产品冲击环境模拟实验的基本参数。 2冲击谱详解 所谓冲击谱，是将冲击源施加于一系列线性、单自由度质量-弹簧系统时，将各单自由度系统的响应运动中的最大响应值，作为对应于系统固有频率的函数而绘制的曲线，即称为冲击谱。由定义可知，冲击谱是单自由度系统受冲击作用后所产生的响应运动在频域中的特性描述。它不同于冲击源的傅里叶频谱，其区别在于：傅里叶频谱仅仅研究冲击源本身在频域中的能量分布属性，只是冲击源函数在频域中的展开，它不涉及任何一个要研究的机械系统的响应。虽然冲击频谱与傅里叶频谱两者都是频率的函数，但有着明显的区别。 换言之，冲击谱是一系列固有频率不同的单自由度线性系统受同一冲击激励响应的总结果。产品受冲击作用，其冲击响应的最大值意味着产品出现最大应力，即试验样品有最大的变形。因此，冲击响应的最大加速度Amax与产品受冲击作用造成的损伤及故障产生的原因直接相关，由此引出了最大冲击响应谱。 3最大冲击响应谱又可以作如下细分 1.正初始冲击响应谱（+I）是指激励脉冲持续时间内，一系列被激励单自由度系统与激励脉冲同方向上出现的最大响应值。Amax(+I)与相应系统的固有频率fn的关系曲线。 2.正残余冲击响应谱（+R）是指激励脉冲持续时间结束后，一系列被激单自由度系统与激励脉冲同方向上出现的最大响应值Amax（+R）与相应系统的固有频率fn的关系曲线。 3.负初始冲击响应谱（-I）是指激励脉冲持续时间内，一系列被激励单自由度系统与激励脉冲反方向上出现的最大值Amax（-I）与相应的系统固有频率fn的关系曲线。 4.负残余冲击响应谱(-R)是指激励脉冲持续时间结束后，一系列被激单自由度系统与激励脉冲反方向上出现的最大值Amax（-R）与相应的系统固有频率fn的关系曲线。 冲击响应谱反映的是环境特性，根据分析冲击响应谱，可以为设计产品的抗冲击能力提供依据。要获取冲击响应谱，首先要采集环境冲击的时域信号，然后再通过软件进行分析，获取冲击响应谱。国内外都有相应的系统可以完成这个工作。比如国内的INTEST（英泰斯特）,提供了冲击加速度时域采集、频域分析、冲击响应谱分析等多种功能，还可以在软件中生成标准脉冲的、归一化后的冲击响应谱，为工程应用提供最直接的解决方案。 4冲击响应谱技术参数 冲击响应谱试验机是用来衡量冲击运动对电工电子产品作用力的大小，考核试品在冲击环境下功能的适应性和结构完好性。 产品特点： 摆锤式结构。 plc控制预设能量自动冲击无二次冲击。 冲击能量无级可调。 计算机测量同时采集时域、频域冲击波形 结合式程序调节器，低频能量调节方便。
功率谱估计介绍 谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。在这里，结合matlab，我做一个粗略介绍。功率谱估计可以分为经典谱估计方法与现代谱估计方法。经典谱估计中最简单的就是周期图法，又分为直接法与间接法。直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积，就得到功率谱的估计；间接法先计算N点样本数据的自相关函数，然后取自相关函数的傅里叶变换，即得到功率谱的估计.都可以编程实现，很简单。在matlab中，周期图法可以用函数periodogram实现。 周期图法估计出的功率谱不够精细，分辨率比较低。因此需要对周期图法进行修正，可以将信号序列x(n)分为n个不相重叠的小段，分别用周期图法进行谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。还可以将信号序列x(n)重叠分段，分别计算功率谱，再计算平均值作为整段数据的功率谱估计。 种称为分段平均周期图法，一般后者比前者效果好。加窗平均周期图法是对分段平均周期图法的改进，即在数据分段后，对每段数据加一个非矩形窗进行预处理，然后在按分段平均周期图法估计功率谱。相对于分段平均周期图法，加窗平均周期图法可以减小频率泄漏，增加频峰的宽度。welch法就是利用改进的平均周期图法估计估计随机信号的功率谱，它采用信号分段重叠，加窗，FFT 等技术来计算功率谱。与周期图法比较，welch法可以改善估计谱曲线的光滑性，大大提高谱估计的分辨率。matlab中，welch法用函数psd实现。调用格式如下： [Pxx,F] = PSD(X,NFFT,Fs,WINDOW,NOVERLAP) X：输入样本数据 NFFT：FFT点数 Fs:采样率 WINDOW：窗类型 NOVERLAP，重叠长度 现代谱估计主要针对经典谱估计分辨率低和方差性不好提出的，可以极大的提高估计的分辨率和平滑性。可以分为参数模型谱估计和非参数模型谱估计。参数模型谱估计有AR模型，MA模型，ARMA模型等；非参数模型谱估计有最小方差法和MUSIC法等。由于涉及的问题太多，这里不再详述，可以参考有关资料。matlab中，现代谱估计的很多方法都可以实现。music方法用pmusic命令实现；pburg函数利用burg法实现功率谱估计；pyulear函数利用yule-walker算法实现功率谱估计等等。 另外，sptool工具箱也具有功率谱估计的功能。窗口化的操作界面很方便，而且有多种方法可以选择 在海杂波抑制的研究中，对海杂波谱分析一定要用到谱估计理论，一定得花时间学好！
常用计算方法 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实现 在应用数学和物理学中，谱密度、功率谱密度和能量谱密度是一个用于信号的通用概念，它表示每赫兹的功率、每赫兹的能量这样的物理量纲。在物理学中，信号通常是波的形式，例如电磁波、随机振动或者声波。当波的频谱密度乘以一个适当的系数后将得到每单位频率波携带的功率，这被称为信号的功率谱密度（power spectral density, PSD）或者谱功率分布（spectral power distribution, SPD）。功率谱密度的单位通常用每赫兹的瓦特数（W/Hz）表示，或者使用波长而不是频率，即每纳米的瓦特数（W/nm）来表示。信号的功率谱密度当且仅当信号是广义的平稳过程的时候才存在。如果信号不是平稳过程，那么自相关函数一定是两个变量的函数，这样就不存在功率谱密度，但是可以使用类似的技术估计时变谱密度。信号功率谱的概念和应用是电子工程的基础，尤其是在电子通信系统中，例如无线电和微波通信、雷达以及相关系统。因此学习如何进行功率谱密度估计十分重要，借助于Matlab工具可以实现各种谱估计方法的模拟仿真并输出结果。下面对周期图法、修正周期图法、最大熵法、Levinson递推法和Burg法的功率谱密度估计方法进行程序设计及仿真并给出仿真结果。 以下程序运行平台：Matlab R2015a（8.5.0.197613） 一、周期图法谱估计程序 1、源程序 Fs=100000; %采样频率100kHz N=1024; %数据长度N=1024 n=0:N-1; t=n/Fs; xn=sin(2000*2*pi*t); %正弦波,f=2000Hz Y=awgn(xn,10); %加入信噪比为10db的高斯白噪声 subplot(2,1,1); plot(n,Y) title('信号') xlabel('时间');ylabel('幅度');
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)
AR模型的谱估计是现代谱估计的主要内容 AR模型的谱估计是现代谱估计的主要内容。 1．AR 模型的Yule—Walker方程和Levinson-Durbin递推算法：在MATLAB中，函数levinson和aryule都采用 Levinson-Durbin递推算法来求解AR模型的参数a1,a2,……,ap及白噪声序列的方差，只是两者的输入参数不同，它们的格式为： A=LEVINSON(R,ORDER) A=ARYULE(x,ORDER) 两函数均为定阶ORDER的求解，但是函数levinson的输入参数要求是序列的自相关函数，而函数aryule的输入参数为采样序列。 下面语句说明函数levinson和函数aryule的功能是相同的： 例子： randn('seed',0) a=[1 0.1 0.2 0.3 0.4 0.5]; x=impz(1,a,20)+randn(20,1)/20; r=xcorr(x,'biased'); r(1:length(x)-1)=[]; A=levinson(r,5) B=aryule(x,5) 2．Burg算法： 格式为：A=ARBURG(x,ORDER); 其中x为有限长序列，参数ORDER用于指定AR 模型的阶数。以上面的例子为例： randn('seed',0) a=[1 0.1 0.2 0.3 0.4 0.5]; x=impz(1,a,20)+randn(20,1)/20; A=arburg(x,5)
3.改进的协方差法： 格式为：A=ARMCOV(x,ORDER); 该函数用来计算有限长序列x(n)的ORDER阶AR 模型的参数。例如：输入下面语句： randn('seed',0) a=[1 0.1 0.2 0.3 0.4 0.5]; x=impz(1,a,20)+randn(20,1)/20; A=armcov(x,5) AR模型阶数P的选择： AR 模型阶数P一般事先是不知道的，需要事先选定一个较大的值，在递推的过程中确定。在使用Levinson—Durbin递推方法时，可以给出由低阶到高阶的每一组参数，且模型的最小预测误差功率Pmin（相当于白噪声序列的方差）是递减的。直观上讲，当预测误差功率P达到指定的希望值时，或是不再发生变化时，这时的阶数即是应选的正确阶数。 因为预测误差功率P是单调下降的，因此，该值降到多少才合适，往往不好选择。比较常见的准则是： 最终预测误差准则：FPE(r)=Pr{[N+(r+1)]/ [N-(r+1)]} 信息论准则：AIC(r)=N*log(Pr)+2*r 上面的N为有限长序列x(n)的长度，当阶数r由1增加时，FPE(r) 和AIC(r)都将在某一r处取得极小值。将此时的r定为最合适的阶数p。 MATLAB中AR模型的谱估计的函数说明： 1． Pyulear函数： 功能：利用Yule--Walker方法进行功率谱估计. 格式： Pxx=Pyulear(x,ORDER,NFFT) [Pxx,W]=Pyulear(x,ORDER,NFFT) [Pxx,W]=Pyulear(x,ORDER,NFFT,Fs) Pyulear(x,ORDER,NFFT,Fs,RANGE,MAGUNITS)
disp(' ') disp(' srs.m ver 2.0 July 3, 2006') disp(' by Tom Irvine Email: tomirvine@https://www.doczj.com/doc/7e13599335.html,') disp(' ') disp(' This program calculates the shock response spectrum') disp(' of an acceleration time history, which is pre-loaded into Matlab.') disp(' The time history must have two columns: time(sec) & acceleration') disp(' ') % clear t; clear y; clear yy; clear n; clear fn; clear a1; clear a2 clear b1; clear b2; clear jnum; clear THM; clear resp; clear x_pos; clear x_neg; % iunit=input(' Enter acceleration unit: 1= G 2= m/sec^2 '); % disp(' ') disp(' Select file input method '); disp(' 1=external ASCII file '); disp(' 2=file preloaded into Matlab '); file_choice = input(''); % if(file_choice==1) [filename, pathname] = uigetfile('*.*'); filename = fullfile(pathname, filename); % fid = fopen(filename,'r'); THM = fscanf(fid,'%g %g',[2 inf]); THM=THM'; else THM = input(' Enter the matrix name: '); end % t=double(THM(:,1));
实验功率谱估计 实验目的: 1、掌握最大熵谱估计的基本原理。 2、了解最终预测误差（FPE）准则。 3、掌握周期图谱估计的基本原理。 4、掌握传统谱估计中直接法与间接法之间的关系。 5、复习快速傅里叶变换与离散傅里叶变换之间关系。 实验内容: 1、设两正弦信号的归一化频率分别为0.175和0.20，用最大熵法编程计算信噪比S/N=30dB、N=32点时该信号的最大熵谱估计结果。 2、用周期图法编程计算上述信号的谱估计结果。 程序示例: 1、最大熵谱估计 clc; N=32; SNR=30; fs=1; t=1:N; t=t/fs; y=sin(2*pi*0.175*t)+sin(2*pi*0.20*t); x = awgn(y,SNR); M=1; P(M)=0; Rx(M)=0; for n=1:N P(M)=P(M)+(abs(x(n)))^2; ef(1,n)=x(n); eb(1,n)=x(n); end P(M)=P(M)/N; Rx(M)=P(M); M=2;
A=0; D=0; for n=M:N A=A+ef(M-1,n)*eb(M-1,n-1); D=D+(abs(ef(M-1,n)))^2+(abs(eb(M-1,n-1)))^2; end xishu=-2*A/D; a(M-1,M-1)=-2*A/D; P(M)=P(M-1)*(1-(abs(xishu))^2); FPE(M-1)=P(M)*(N+M)/(N-M); TH=FPE(M-1); for n=M:N ef(M,n)=ef(M-1,n)+xishu*eb(M-1,n-1); eb(M,n)=eb(M-1,n-1)+xishu*ef(M-1,n); end M=M+1; A=0; D=0; for n=M:N A=A+ef(M-1,n)*eb(M-1,n-1); D=D+(abs(ef(M-1,n)))^2+(abs(eb(M-1,n-1)))^2; end xishu=-2*A/D; a(M-1,M-1)=-2*A/D; P(M)=P(M-1)*(1-(abs(xishu))^2); FPE(M-1)=P(M)*(N+M)/(N-M); for m=1:M-2 a(M-1,m)=a(M-2,m)+xishu*a(M-2,M-1-m); end while FPE(M-1)
【关键字】技术 冲击响应谱校准技术的研究 厉巍陈永久朱永晓 （贵州航天计量测试技术研究所，贵州贵阳550009） 摘要：冲击响应谱试验已经成为大多数航天产品必做的力学环境试验项目之一，传统的冲击试验缺乏对冲击环境模拟的真实性，本文介绍了冲击响应谱的原理和冲击响应谱试验设备；用labVIEW为平台，编写了冲击响应谱校准软件，为冲击响应谱试验机的校准与数据分析提供了通用性较好的校准分析方法，并基于PXI系统设计了冲击响应谱校准装置。 关键词：航天产品LabVIEW 冲击响应谱校准PXI系统 0引言 冲击响应谱试验机是用于完成冲击响应谱试验的环境试验设备，冲击响应谱是对产品实施抗冲击设计的分析基础，也是生产中冲击环境模拟试验的基本参数，在航空、航天重点型号科研生产及有关重大科技专项中，冲击响应谱试验已经成为必做的环境试验之一。产品在实际应用过程中受力情况复杂，其中，冲击激励会使设备激起强迫振动和固有频率响应，使产品性能和结构强度受到不同程度的损害甚至失效。航空、航天、电子等行业产品在生产、运输等过程中存在着各种冲击，而这对产品的质量和可靠性有着很大的负面影响。为了解决这一问题，在此基础上产生并发展起了冲击试验。近年来，随着对环境试验的认识不断提高，对冲击环境的模拟也提出了更高的要求，冲击响应谱试验也来越被关注。 1 冲击响应谱原理 冲击信号与一般的振动信号在许多方面具有不同的特性，工程中研究冲击信号的目的并不是研究冲击波形本身，而是更加注重冲击作用于系统的效果，或者说是研究冲击运动对系统的损伤势。不论用冲击的时间历程还是用频谱都难以描述冲击的损伤势，因此必须使用能够衡量冲击效果的冲击响应谱。 冲击响应谱系指一单自由度质量弹簧阻尼系统，当公共基础受到冲击激励时产生的响应峰值作为单自由度系统固有频率的函数绘出的图，其物理模型如图1所示。 图1 冲击响应谱的物理模型 数学模型可归结为如下微分方程的解： 式中，； ； 2 冲击响应谱试验设备
matlab实现经典功率谱估计 fft做出来是频谱，psd做出来是功率谱；功率谱丢失了频谱的相位信息；频谱不同的信号其功率谱是可能相同的；功率谱是幅度取模后平方，结果是个实数 matlab中自功率谱密度直接用psd函数就可以求，按照matlab的说法，psd能实现Welch法估计，即相当于用改进的平均周期图法来求取随机信号的功率谱密度估计。psd求出的结果应该更光滑吧。 1、直接法： 直接法又称周期图法，它是把随机序列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)); 2、间接法： 间接法先由序列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;
冲击响应谱(SRS)是一个瞬态加速度脉冲可能对结构造成破坏的图示。它绘制了一组单自由度(SDOF)弹簧的峰值加速度响应，就像在刚性无质量的基础上一样,质量阻尼器系统都经历相同的基本激励。每个SDOF系统具有不同的固有频率;它们都有相同的粘滞阻尼因子。频谱的结果是在固有频率(水平方向)上绘制峰值加速度(垂直)得出的。一个SRS是由一个冲击波产生，使用以下过程: 指定SRS的阻尼比(5%是最常见的)、使用数字滤波器模拟频率单自由度、fn和阻尼ξ。应用瞬态作为输入，计算响应加速度波形。保留在脉冲持续时间和之后的峰值正负响应。选择其中一个极值，并将其绘制成fn的频谱振幅。对每个(对数间隔)fn重复这些步骤。 由此产生的峰值加速度与弹簧-质量阻尼系统固有频率的曲线称为冲击响应谱，简称SRS。一个SDOF机械系统由以下组件组成: ①质量,米 ②弹簧,K ③阻尼器,C Fn,固有频率和临界阻尼因子,ξ,描述一个应用系统,可以从上面的参数计算。对于小于或等于0.05的小阻尼比，频率响应的峰值发生在fn的邻近区域，其中
Q为质量因子，等于1/(2ξ)。 任何瞬态波形都可以作为SRS呈现，但这种关系不是唯一的;许多不同的瞬态波形可以产生相同的SRS。SRS不包含所有关于瞬态波形的信息，因为它只跟踪峰值瞬时加速度。 不同的阻尼比为相同的冲击波形产生不同的SRS。零阻尼会产生最大的响应，而高阻尼则会产生较平的SRS。阻尼比与质量因子Q有关，在正弦振动的情况下也可以被认为是传递率。阻尼比为5%（ξ=0.05）时，Q值为10。如果没有指定阻尼因子(或Q)，则SRS图是不完整的。 ★SRS箱的频率间隔 一个SRS由多个在对数频率范围内均匀分布的箱组成。频率分布可以由两个数字来定义:一个参考频率和期望的分数倍频间隔，如1/1、1/3或1/6。(倍频程是频率的两倍)例如，250hz和500hz的频率相差一个倍频程， 1 kHz和2 kHz的频率也是一样。 比例带宽显示对于分析各种自然系统，如人类对噪声和振动的反应，是非常有用的。许多机械系统表现出的特征非常适合以比例带宽分析。 为了获得更好的频率分辨率，频率范围可以以倍频程的一部分划分比例间隔。例如，有1/3倍频间隔，每个倍频程有3个SDOF滤波器。一般来说，对于1/N个分数倍频程，每个倍频程有N个带通滤波器。这里1/N称为分数倍
《计算方法》实验报告 指导教师： 学院： 班级： 团队成员：
一、题目 例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)掌握三种算法的概念、应用及特点； (2)了解谱估计在信号分析中的作用； (3) 能够利用burg 法对信号作谱估计，对信号的特点加以分析。 二;实验内容： （1）简单说明三种方法的原理。 （2）用三种方法编写程序，在matlab 中实现。 （3）将计算结果表示成图形的形式，给出三种情况的功率谱图。 （4）比较三种方法的特性。 （5）写出自己的心得体会。 三：实验原理： 1.周期图法： 周期图法又称直接法。它是从随机信号x(n)中截取N 长的一段，把它视为能量有限x(n)真实功率谱)(jw x e S 的估计)(jw x e S 的抽样. 认为随机序列是广义平稳且各态遍历的，可以用其一个样本x(n)中的一段)(n x N 来估计该随机序列的功率谱。这当然必然带来误差。由于对)(n x N 采用DFT ，就默认)(n x N 在时域是周期的，以及)(k x N 在频域是周期的。这种方法把随机序列样本x(n)看成是截得一段)(n x N 的周期延拓，这也就是周期图法这个名字的来历。
2.相关法（间接法): 这种方法以相关函数为媒介来计算功率谱，所以又叫间接法。这种方法的具体步骤是： 第一步：从无限长随机序列x(n)中截取长度N 的有限长序列列 )(n x N 第二步：由N 长序列)(n x N 求（2M-1）点的自相关函数)(m R x ∧ 序列。 )()(1)(1 m n x n x N m R N n N N x += ∑-=∧ (2-1) 这里，m=-(M-1)…,-1,0,1…,M-1，M N ，)(m R x 是双边序列，但是由自相关函数的偶对称性式，只要求出m=0,。。。,M-1的傅里叶变换，另一半也就知道了。 第三步：由相关函数的傅式变换求功率谱。即 jwm M M m X jw x e m R e S ----=∧∧ ∑= )()(1) 1( 以上过程中经历了两次截断，一次是将x(n)截成N 长，称为加数据窗，一次是将x(n)截成（2M-1）长，称为加延迟窗。因此所得的功率谱仅是近似值，也叫谱估计，式中的)(jw x e S 代表估值。一般取M<
相关主题
文本预览