《数据采集与信号处理作业基于FFT的功率谱分析程序设计与应用》由会员分享,可在线阅读,更多相关《数据采集与信号处理作业基于FFT的功率谱分析程序设计与应用(15页珍藏版)》请在金锄头文库上搜索。
1、数据采集与信号处理作业一、基本内容: 基于FFT的功率谱分析程序设计与应用1基本要求 1) 对一个人为产生的信号进行采用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计算。处理结果为采用窗口显示时域波形和频域波形。基于MATLAB6.5的源程序如下:clf;clc;Fs=300;%采样频率n=0:1:300;N=128;Sf=20; %信号频率xn=120*cos(2*pi*Sf*n/Fs);%产生波形序列cxn=xcorr(xn,unbiased); %计算序列的自相关函数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)信号经过均值化处理或不经过均值化处理的结果比较%含直流分量而未均值化的信号clc;clf;Fs=300; %采样频率n=0:1:300;N=128;Sf=20; %信号频率xn=120*cos(2*pi*Sf*n/Fs)+60;%产生波形序列window=boxcar(length(xn); %矩形窗nfft=512;%采样点数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;%产生波形序列z=mean(xn);h=xn-z;window=boxcar(length(xn); %矩形窗nfft=512;%采样点数Pxx,f=periodogram(h
6、,window,nfft,Fs); %直接法subplot(2,1,2);plot(f,Pxx);title(含有直流分量的余弦曲线均值化后的功率谱波形图);xlabel(Frequency(Hz);ylabel(Magnitude);余弦函数信号x的直流分量即x的均值,它的存在将导致余弦函数在f0 (即图形的波峰处的频率)处的功率降低,因此,必须对带有直流分量的信号进行均值化处理,以消除直流分量的干扰。由图3可以看出,带有直流分量的余弦函数不经均值化处理时的功率谱在低频段含有一个很大的功率分量,导致余弦函数在f0处的功率(或主频率分量)降低。而经均值化处理后,余弦函数在f0处的功率明显增加,
7、这一点可通过比较上述两图中f0处的波峰值的大小得以证明。2)采用不同窗函数时的谱结果(1)余弦函数加矩形窗的功率谱图显示如下:(2)余弦函数加汉宁的功率谱图显示如下:(3)余弦函数加汉明窗的功率谱图显示如下:(4)余弦函数加布拉克曼明窗的功率谱图显示如下:MATLAB程序如下:%矩形框,汉宁框,汉明框,布拉克曼框函数处理结果;clc;Fs=300;%采样频率Sf=20;t=0:1/Fs:1;xn=120*cos(2*pi*Sf*t);%产生信号序列nfft=1024;window1=boxcar(length(xn);%矩形窗函数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、)典型函数的频谱(1)矩形窗函数的时域频域波形图:(2)汉宁窗函数的时域,频域波形图:(3)直线的频域波形图:(4)跃函数的时域,频域波形图(5)函数的时域,频域波形图(6)方波函数的时域,频域波形图(7)三角波函数的时域,频域波形图相关程序如下:%典型波形产生(矩形窗函数, 汉宁窗函数,直线,阶跃函数,函数,方波,三角波等)%矩形窗函数t=0:0.001:0.2;N=256;FS=300;w=boxcar(N); %产生信号subplot(211);plot(w);title(矩形窗函数的时域波形图);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(矩形窗函数频域波形图);grid;xlabel(Frequency(Hz);ylabel(Magnitude);%汉宁窗函数的频域波形图t=0:0.001:0.2;N=256;FS=300;w=hanning(N); %产生信号subplot(211);plot(w);title(汉宁窗函数的时域波形图);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(汉宁窗函数频域波形图);grid;xlabel(Frequency(Hz);ylabel(Magnitude);%直线的频域波形图t=0:0.001:0.2;N=256;FS=300;w=1; %产生信号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);%阶跃函数的频域波形图clc;clf;t=0:0.001:0.2;N=256;FS=300;w=ones(1,N); %产生信号subplot(211);plot(w);title(阶跃函数的时域波形图);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;xlabel(Frequency(Hz);ylabel(Magnitude);%函数的频域波形图clc;clf;t=0:0.001:0.2;N=256;FS=300;w=zeros(1,N);w(1)=1; %产生信号subplot