对于频谱分析能量守恒的验证Word文档下载推荐.docx

对于频谱分析能量守恒的验证Word文档下载推荐.docx

《对于频谱分析能量守恒的验证Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《对于频谱分析能量守恒的验证Word文档下载推荐.docx（19页珍藏版）》请在冰豆网上搜索。

xinhaohz2=50*df;

f1=3;

f2=4;

x=f1*sin（2*pi*xinhaohz1*t）+f2*sin（2*pi*xinhaohz2*t）;

m=8;

N=256;

nxd=bin2dec（fliplr（dec2bin（[1:

N]-1,m）））+1;

y=x（nxd）;

formm=1:

m

Nz=2^mm;

u=1;

WN=exp（-i*2*pi/Nz）;

forj=1:

Nz/2

fork=j:

Nz:

N

kp=k+Nz/2;

t=y（kp）*u;

y（kp）=y（k）-t;

y（k）=y（k）+t;

end

u=u*WN;

end

%自己编的FFT跟直接调用的函数运算以后的结果进行对比

y1=fft（x,256）;

Pyy=（abs（y）/nfft）.^2;

f=df*（0:

127）;

figure

（1）;

subplot（1,2,1）;

plot（f,Pyy（1:

128））

axis（[0,600,0,4.5]）

title（'

.mbyhand'

）

xlabel（'

Frequency（Hz）'

Pyy=（abs（y1）/nfft）.^2;

subplot（1,2,2）;

.mbylibraryfunction'

%fftHvsL.m程序结束，MATLAB库函数fft（）和手工函数对信号的频谱分析结束

执行结果：

图表1MATLAB库函数fft（）和手工函数对的频谱分析，频谱数据见文档“频谱数据.txt”

结论：

手工fft函数是正确的，和MATLABfft库函数一致

B.环境：

ampDH.m程序：

%ampDH.m程序开始，MATLAB库函数实现FFT的振幅检测

tfinal=dt*nfft;

x=3*sin（2*pi*xinhaohz1*t）+4*sin（2*pi*xinhaohz2*t）;

plot（x（1:

50））

y=fft（x,nfft）;

E=sum（Pyy）

figure

（2）;

Powerspectraldensity'

%ampDL.m程序结束，MATLAB库函数实现FFT的振幅检测结束

%ampDH.m程序开始，MATLAB中用手工函数实现FFT的振幅检测

%y=fft（x,nfft）;

bit_rev=[

0,128,64,192,32,160,96,224,16,144,80,208,48,176,112,240,8,136,72,200,40,168,104,232,24,152,88,216,56,184,120,248,4,132,68,196,36,164,100,228,20,148,84,212,52,180,116,244,12,140,76,204,44,172,108,236,28,156,92,220,60,188,124,252,2,130,66,194,34,162,98,226,18,146,82,210,50,178,114,242,10,138,74,202,42,170,106,234,26,154,90,218,58,186,122,250,6,134,70,198,38,166,102,230,22,150,86,214,54,182,118,246,14,142,78,206,46,174,110,238,30,158,94,222,62,190,126,254

];

foritemp=1:

128

bit_revTemp=bit_rev（itemp）;

y（itemp）=x（bit_revTemp+1）;

y（257-itemp）=x（256-bit_revTemp）;

end;

U=complex（0）;

W=complex（0）;

T=complex（0）;

LE=1;

LE1=1;

I=1;

J=1;

IP=1;

L=1;

while（L<

=8）

LE=LE*2;

LE1=LE/2;

U=complex（1,0）;

W=complex（cos（pi/（1.0*LE1））,-1.0*sin（pi/（1.0*LE1）））;

J=1;

while（J<

=LE1）

I=J;

while（I<

=256）

IP=I+LE1;

T=complex（real（y（IP））*real（U）-imag（y（IP））*imag（U）,real（y（IP））*imag（U）+imag（y（IP））*real（U））;

y（IP）=complex（real（y（I））-real（T）,imag（y（I））-imag（T））;

y（I）=complex（real（y（I））+real（T）,imag（y（I））+imag（T））;

I=I+LE;

end;

tempfloat=real（U）;

U=complex（real（U）*real（W）-imag（U）*imag（W）,tempfloat*imag（W）+imag（U）*real（W））;

J=J+1;

L=L+1;

%ampDH.m程序结束，MATLAB中用手工函数实现FFT的振幅检测结束

程序执行结果：

图表2MATLAB库函数实现FFT的频谱图E=12.5000

图表3MATLAB中用手工函数实现FFT的频谱图E=12.5000

手工fft函数和MATLAB库函数fft（）实现的振动幅值和频率检测是一致的

C.环境：

DiDH.m程序：

%DiDH.m程序开始，用手工函数实现的通过频域积分来计算位移均方值

f1=50;

f2=40;

%加速度信号是x，位移信号是d——x=（d（d））^2/dt^2；

对应一个对t的二次积分

d=（-f1/（（2*xinhaohz1*pi）^2））*sin（2*pi*xinhaohz1*t）+（-f2/（（2*pi*xinhaohz2）^2））*sin（2*pi*xinhaohz2*t）;

%djunfangzhi为位移信号均方值；

dfengzhi为位移峰值，简谐信号的均方值和峰值关系是：

峰值=sqrt（均方值*2）

djunfangzhi=0;

fori=1:

nfft

djunfangzhi=djunfangzhi+d（i）^2;

djunfangzhi=djunfangzhi/nfft

dfengzhi=sqrt（djunfangzhi）/0.707

%海明窗系数矩阵

hamm=[

0.0800,0.0801,0.0806,0.0813,0.0822,0.0835,0.0850,0.0868,0.0889,0.0913,0.0939,0.0968,0.1000,0.1034,0.1071,0.1111,0.1153,0.1198,0.1245,0.1295,0.1347,0.1402,0.1459,0.1519,0.1581,0.1645,0.1712,0.1781,0.1852,0.1925,0.2001,0.2078,0.2157,0.2239,0.2322,0.2407,0.2494,0.2583,0.2673,0.2765,0.2859,0.2954,0.3051,0.3149,0.3249,0.3350,0.3452,0.3555,0.3659,0.3765,0.3871,0.3979,0.4087,0.4196,0.4305,0.4416,0.4527,0.4638,0.4750,0.4863,0.4976,0.5089,0.5202,0.5315,0.5428,0.5542,0.5655,0.5768,0.5881,0.5993,0.6106,0.6217,0.6329,0.6439,0.6549,0.6659,0.6767,0.6875,0.6982,0.7088,0.7193,0.7297,0.7400,0.7501,0.7601,0.7700,0.7797,0.7893,0.7988,0.8081,0.8172,0.8262,0.8350,0.8436,0.8520,0.8602,0.8683,0.8761,0.8837,0.8912,0.8984,0.9054,0.9121,0.9187,0.9250,0.9311,0.9369,0.9426,0.9479,0.9530,0.9579,0.9625,0.9669,0.9710,0.9748,0.9784,0.9817,0.9847,0.9875,0.9899,0.9922,0.9941,0.9958,0.9972,0.9983,0.9991,0.9997,1.0000

%FFT计算前将时域加速度序列进行逆序输入的矩阵

0,128,64,192,32,160,96,224,16,144,80,208,48,176,112,240,8,136,72,20