數(shù)據(jù)采集與信號處理作業(yè)基于FFT的功率譜分析程序設計與應用

1、數(shù)據(jù)采集與信號處理作業(yè)一、基本內(nèi)容： 基于fft的功率譜分析程序設計與應用1基本要求 1） 對一個人為產(chǎn)生的信號進行采用fft變換方法進行功率譜分析。 已知信號x(n)=120.0*cos(2*3.14*sf*n/fs) 式中： n=0,1,2 n-1 sf-信號頻率 fs-采樣頻率 其fft變換結果x(k)可用下面提供的fft子程序求出，計算功率譜的公式為： w（k）=2(xr(k)2 +xi(k)2)/n 式中： k=0,1,2 n/2-1 xr(k)- x(k)的實部 xi(k)- x(k)的虛部 請用vb,vc或c+builder編譯器編程，或采用matlab計算，或采用高級語言調用m

2、atlab計算。處理結果為采用窗口顯示時域波形和頻域波形?；趍atlab6.5的源程序如下：clf;clc;fs=300;%采樣頻率n=0:1:300;n=128;sf=20; %信號頻率xn=120*cos(2*pi*sf*n/fs);%產(chǎn)生波形序列cxn=xcorr(xn,unbiased); %計算序列的自相關函數(shù)figure(1);plot(n,xn);grid on;title(余弦曲線的時頻圖);nfft=1024;cxk=fft(cxn,nfft);pxx=abs(cxk); %求序列的psd;index=0:round(nfft/2-1);k=index*fs/nfft;pl

3、ot_pxx=pxx(index+1);figure(2);plot(k,plot_pxx);title(余弦曲線的功率譜曲線圖);xlabel(frequency(hz);ylabel(magnitude);grid on; 2）對實驗所采集的轉子振動信號進行頻譜分析sf=1000; %采樣頻率1024/508hzfid = fopen(d:fanbo_45hz_1024_1000hzsanjiao_45hz_1024_1000hz);%轉子信號a,n= fscanf(fid,%f); fclose(fid);y=fft(a,n);%fft運算pyy =sqrt(y.*conj(y)*2.0

4、/n;%取功率普密度f=(0:length(pyy)-1)*sf/length(pyy);lpyy=20*log10(pyy);plot(f(1:n/2),pyy(1:n/2),black);%輸出fs/2點幅頻譜圖grid;2.討論1）信號經(jīng)過均值化處理或不經(jīng)過均值化處理的結果比較%含直流分量而未均值化的信號clc;clf;fs=300; %采樣頻率n=0:1:300;n=128;sf=20; %信號頻率xn=120*cos(2*pi*sf*n/fs)+60;%產(chǎn)生波形序列window=boxcar(length(xn); %矩形窗nfft=512;%采樣點數(shù)pxx,f=periodogra

5、m(xn,window,nfft,fs); %直接法subplot(2,1,1);plot(f,pxx);title(含有直流分量的余弦曲線未均值化的功率譜波形圖);xlabel(frequency(hz);ylabel(magnitude);%含直流分量且均值化的信號clc;clf;fs=300; %采樣頻率n=0:1:300;n=128;sf=20; %信號頻率xn=120*cos(2*pi*sf*n/fs)+60;%產(chǎn)生波形序列z=mean(xn);h=xn-z;window=boxcar(length(xn); %矩形窗nfft=512;%采樣點數(shù)pxx,f=periodogram(h

6、,window,nfft,fs); %直接法subplot(2,1,2);plot(f,pxx);title(含有直流分量的余弦曲線均值化后的功率譜波形圖);xlabel(frequency(hz);ylabel(magnitude);余弦函數(shù)信號x的直流分量即x的均值，它的存在將導致余弦函數(shù)在f0 (即圖形的波峰處的頻率)處的功率降低，因此，必須對帶有直流分量的信號進行均值化處理，以消除直流分量的干擾。由圖3可以看出，帶有直流分量的余弦函數(shù)不經(jīng)均值化處理時的功率譜在低頻段含有一個很大的功率分量，導致余弦函數(shù)在f0處的功率（或主頻率分量）降低。而經(jīng)均值化處理后，余弦函數(shù)在f0處的功率明顯增加，

7、這一點可通過比較上述兩圖中f0處的波峰值的大小得以證明。2）采用不同窗函數(shù)時的譜結果（1）余弦函數(shù)加矩形窗的功率譜圖顯示如下：（2）余弦函數(shù)加漢寧的功率譜圖顯示如下：（3）余弦函數(shù)加漢明窗的功率譜圖顯示如下：（4）余弦函數(shù)加布拉克曼明窗的功率譜圖顯示如下：matlab程序如下：%矩形框，漢寧框，漢明框，布拉克曼框函數(shù)處理結果；clc;fs=300;%采樣頻率sf=20;t=0:1/fs:1;xn=120*cos(2*pi*sf*t);%產(chǎn)生信號序列nfft=1024;window1=boxcar(length(xn);%矩形窗函數(shù)window2=hanning(length(xn);%漢寧窗w

8、indow3=hamming(length(xn);%漢明窗window4=blackman(length(xn);%布拉克曼窗noverlap=0;p=0.9;pxx1,f1=psd(xn,nfft,fs,window1,noverlap,p);%計算序列的psd;index=0:round(nfft/2-1);k=index*fs/nfft;plot_pxx1=pxx1(index+1);figure(1)plot(k,plot_pxx1);gridtitle(矩形窗);pxx2,f2=psd(xn,nfft,fs,window2,noverlap,p);%計算序列的psd;plot_px

9、x2=pxx2(index+1);figure(2)plot(k,plot_pxx2);gridtitle(漢寧窗);pxx3,f3=psd(xn,nfft,fs,window3,noverlap,p);%計算序列的psd;plot_pxx3=pxx3(index+1);figure(3)plot(k,plot_pxx3);gridtitle(漢明窗);pxx4,f4=psd(xn,nfft,fs,window4,noverlap,p);%計算序列的psd;plot_pxx4=pxx4(index+1);figure(4)plot(k,plot_pxx4);gridtitle(布拉克曼窗);3

10、)典型函數(shù)的頻譜（1）矩形窗函數(shù)的時域頻域波形圖：（2）漢寧窗函數(shù)的時域，頻域波形圖：（3）直線的頻域波形圖：（4）躍函數(shù)的時域，頻域波形圖（5）函數(shù)的時域，頻域波形圖（6）方波函數(shù)的時域，頻域波形圖（7）三角波函數(shù)的時域，頻域波形圖相關程序如下：%典型波形產(chǎn)生(矩形窗函數(shù), 漢寧窗函數(shù)，直線，階躍函數(shù)，函數(shù)，方波，三角波等)%矩形窗函數(shù)t=0:0.001:0.2;n=256;fs=300;w=boxcar(n); %產(chǎn)生信號subplot(211);plot(w);title(矩形窗函數(shù)的時域波形圖);axis(0,260,0,2);grid on;y=fft(w,n); %fft運算mag

11、=abs(y); %取幅值f=(0:length(y)-1)*fs/length(y);subplot(212);plot(f(1:n/2),mag(1:n/2); %輸出fs/2點幅頻譜圖title(矩形窗函數(shù)頻域波形圖);grid;xlabel(frequency(hz);ylabel(magnitude);%漢寧窗函數(shù)的頻域波形圖t=0:0.001:0.2;n=256;fs=300;w=hanning(n); %產(chǎn)生信號subplot(211);plot(w);title(漢寧窗函數(shù)的時域波形圖);grid on;y=fft(w,n); %fft運算mag=abs(y); %取幅值f=(

12、0:length(y)-1)*fs/length(y);subplot(212);plot(f(1:n/2),mag(1:n/2); %輸出fs/2點幅頻譜圖title(漢寧窗函數(shù)頻域波形圖);grid;xlabel(frequency(hz);ylabel(magnitude);%直線的頻域波形圖t=0:0.001:0.2;n=256;fs=300;w=1; %產(chǎn)生信號y=fft(w,n); %fft運算mag=abs(y); %取幅值f=(0:length(y)-1)*fs/length(y);plot(f(1:n/2),mag(1:n/2); %輸出fs/2點幅頻譜圖title(直線頻域

13、波形圖);grid;xlabel(frequency(hz);ylabel(magnitude);%階躍函數(shù)的頻域波形圖clc;clf;t=0:0.001:0.2;n=256;fs=300;w=ones(1,n); %產(chǎn)生信號subplot(211);plot(w);title(階躍函數(shù)的時域波形圖);grid;y=fft(w,n); %fft運算mag=abs(y); %取幅值f=(0:length(y)-1)*fs/length(y);subplot(212);plot(f(1:n/2),mag(1:n/2); %輸出fs/2點幅頻譜圖title(階躍函數(shù)的頻域波形圖);grid;xlab

14、el(frequency(hz);ylabel(magnitude);%函數(shù)的頻域波形圖clc;clf;t=0:0.001:0.2;n=256;fs=300;w=zeros(1,n);w(1)=1; %產(chǎn)生信號subplot(211);plot(w);grid;title(函數(shù)的時域波形圖);y=fft(w,n); %fft運算mag=abs(y); %取幅值f=(0:length(y)-1)*fs/length(y);subplot(212);plot(f(1:n/2),mag(1:n/2); %輸出fs/2點幅頻譜圖title(函數(shù)的頻域波形圖);grid;xlabel(frequency

15、(hz);ylabel(magnitude);%方波函數(shù)的頻域波形圖clc;clf;t=0:0.001:0.2;n=256;fs=300;w=square(2*pi*50*t); %產(chǎn)生信號subplot(211);plot(t,w); title(方波的時域波形圖);axis(0,0.2,-0.2,1.2);grid; y=fft(w,n); %fft運算mag=abs(y); %取幅值f=(0:length(y)-1)*fs/length(y);subplot(212);plot(f(1:n/2),mag(1:n/2); %輸出fs/2點幅頻譜圖title(方波的頻域波形圖);grid;x

16、label(frequency(hz);ylabel(magnitude);%三角波的頻域波形圖t=0:0.001:0.2;n=256;fs=300;w=sawtooth(2*pi*50*t,0.5); subplot(211);plot(t,w); grid;title(三角波的時域波形圖);%產(chǎn)生信號y=fft(w,n); %fft運算mag=abs(y); %取幅值f=(0:length(y)-1)*fs/length(y);subplot(212);plot(f(1:n/2),mag(1:n/2); %輸出fs/2點幅頻譜圖title(三角波的頻域波形圖);grid;xlabel(fr

17、equency(hz);ylabel(magnitude);4）整周期和非整周期采樣時兩者的比較%漢明窗周期采樣;fs=300; %采樣頻率n=0:1:300;n=256;sf=20; %信號頻率xn=120*cos(2*pi*sf*n/fs);nfft=1024;%采樣點數(shù)window3=hamming(length(xn);%漢明框pxx3,f3=periodogram(xn,window3,nfft,fs); %直接法subplot(211);plot(f3,pxx3);grid;title(余弦信號周期采樣);xlabel(frequency(hz);ylabel(magnitude)

18、;%漢寧窗非周期采樣;t=1.7*n; %設置非周期采樣x1=0.5*(1-cos(2*pi*t/n); %定義窗函數(shù)z=xn.*x1;y=fft(z,n); %fft運算mag=abs(y); %取幅值f=(0:length(y)-1)*fs/length(y);subplot(212);plot(f(1:n/2),mag(1:n/2); %輸出fs/2點幅頻譜圖xlabel(frequency(hz);ylabel(magnitude);title(余弦信號非周期采樣);grid;二、設計計算emd分析程序，針對一個模擬調幅信號+加隨機噪聲，選擇需要的imf單個分量ci(t)進行包絡譜分析已知模擬信號為：x=80.*cos(2*pi*f1.*t).*(1+0.8.*cos(2*pi*f