数据采集与信号处理.doc

数据信息采集与处理基本内容：基于FFT的功率谱分析程序设计与应用1. 基本要求1) 对一个人为产生的信号进行采用FFT变换方法进行功率谱分析。已知信号x(n)=80.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计算，或采用高级语言调用MATLAB计算。处理结果为采用窗口显示时域波形和频域波形。此信号的时域谱、频域谱、功率谱如下面图1图3所示：图1图2图3其MATLAB代码为：FS=200;SF=10;N=1024;n=0:N-1;t=n/FS;x=80.0*cos(2*3.14*SF*t);figure;plot(t，x);xlabel(t);ylabel(y);title(x=80.0*cos(2*3.14*SF*t)时域波形);grid;y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换figure;plot(f(1:N/2),mag(1:N/2);%做频谱图xlabel(频率(Hz);ylabel(幅值);title(x=80.0*cos(2*3.14*SF*t)幅频谱图N=1024);grid;Py =2*(y.*conj(y)/N; %计算功率谱密度Pyfigure;plot(f(1:N/2),Py(1:N/2);xlabel(频率(Hz);ylabel(功率谱密度);title(x=80.0*cos(2*3.14*sf*t)功率谱密度);grid;2) 对实验所采集的转子振动信号进行频谱分析图4其MATLAB代码：SF=1000; fid = fopen(D:数据采集与信号处理作业信号处理大作业sanjiao_45HZ_1024_1000HZ);%转子信号a,N= fscanf(fid,%f); fclose(fid);y=fft(a,N);%FFT运算Pyy =sqrt(y.*conj(y)*2.0/N;%取功率普密度f=(0:length(Pyy)-1)*SF/length(Pyy);LPyy=20*log10(Pyy);plot(f(1:N/2),Pyy(1:N/2);%输出FS/2点幅频谱图xlabel(频率(Hz);ylabel(幅值);title(转子振动信号频谱图);grid;2. 讨论1) 信号经过均值化处理或不经过均值化处理的结果比较：图5图6通过以上两个图（图5、图6）的分析，我们可以看出均值化处理后的频谱的低频段消失，这就去去除了常规的干扰频谱，如环境噪声等，对我们进行频谱分析有很大作用。其MATLAB代码为：Fs=200; %采样频率n=0:1:200;N=1024;SF=10; %信号频率xn=80.0*cos(2*3.14*SF*n/Fs)+70;%产生波形序列window=boxcar(length(xn); %矩形窗nfft=512;%采样点数Pxx,f=periodogram(xn,window,nfft,Fs); %直接法plot(f,Pxx);xlabel(频率(Hz);ylabel(功率谱密度);title(含有直流分量的余弦曲线未均值化的功率谱波形图);grid;2) 采用不同窗函数时的谱结果（矩形窗函数, 汉宁窗函数，汉明窗）如图712所示。图7图8图9图10图11图12其MATLAB代码为：y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换w_han=(hanning(N);y1=x.*w_han;figure;plot(t,y1);xlabel(t);ylabel(y);title(汉宁窗时域波形);grid;y2=mag.*w_han;figure;plot(f(1:N/2),y2(1:N/2);xlabel(频率(Hz);ylabel(幅值);title(汉宁窗频域特性);grid;w_rect=(rectwin(N);y3=x.*w_rect;figure;plot(t,y3);xlabel(t);ylabel(y);title(矩形窗时域波形);grid;y4=mag.*w_rect;figure;plot(f(1:N/2),y4(1:N/2);xlabel(频率(Hz);ylabel(幅值);title(矩形窗频域特性);grid;w_ham=(hamming(N);y5=x.*w_ham;figure;plot(t,y5);xlabel(t);ylabel(y);title(汉明窗时域波形);grid;y6=mag.*w_ham;figure;plot(f(1:N/2),y6(1:N/2);xlabel(频率(Hz);ylabel(幅值);title(汉明窗频域特性);grid;3) 典型函数的频谱（矩形窗函数, 汉宁窗函数，直线，阶跃函数，函数，方波，三角波等），如图1318所示。图13图 14图 15图 16图 17图18此部分MATLAB代码如下：t=0:0.001:0.2;N=256;FS=300;w=boxcar(N); %产生信号figure;plot(w);title(矩形窗函数的时域波形图);axis(0,260,0,2);grid on;y=fft(w,N); %FFT运算mag=abs(y);%取幅值f=(0:length(y)-1)*FS/length(y);figure;plot(f(1:N/2),mag(1:N/2); %输出FS/2点幅频谱图title(矩形窗函数频域波形图);grid;xlabel(频率);ylabel(幅值);t=0:0.001:0.2;N=256;FS=300;w=hanning(N); %产生信号figure;plot(w);title(汉宁窗函数的时域波形图);grid on;y=fft(w,N); %FFT运算mag=abs(y); %取幅值f=(0:length(y)-1)*FS/length(y);figure;plot(f(1:N/2),mag(1:N/2); %输出FS/2点幅频谱图title(汉宁窗函数频域波形图);grid on;xlabel(频率);ylabel(幅值);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);figure;plot(f(1:N/2),mag(1:N/2); %输出FS/2点幅频谱图title(直线频域波形图);grid on;xlabel(Frequency(Hz);ylabel(Magnitude);%阶跃函数的频域波图clc;clf;t=0:0.001:0.2;N=256;FS=300;w=ones(1,N); %产生信号figure;plot(w);title(阶跃函数的时域波形图);grid on;y=fft(w,N); %FFT运算mag=abs(y);%取幅值f=(0:length(y)-1)*FS/length(y);figure;plot(f(1:N/2),mag(1:N/2); %输出FS/2点幅频谱图title(阶跃函数的频域波形图);grid on;xlabel(频率);ylabel(幅值);t=0:0.001:0.2;N=256;FS=300;w=zeros(1,N);w(1)=1; %产生信号figure;plot(w);grid on;title(函数的时域波形图);y=fft(w,N);%FFT运算mag=abs(y);%取幅值f=(0:length(y)-1)*FS/length(y);figure;plot(f(1:N/2),mag(1:N/2); %输出FS/2点幅频谱图title(函数的频域波形图);grid on;xlabel(频率);ylabel(幅值);t=0:0.001:0.2;N=256;FS=300;w=square(2*pi*50*t); %产生信号figure;plot(t,w); title(方波的时域波形图);axis(0,0.2,-0.2,1.2);grid on; y=fft(w,N); %FFT运算mag=abs(y); %取幅值f=(0:length(y)-1)*FS/length(y);figure;plot(f(1:N/2),mag(1:N/2); %输出FS/2点幅频谱图title(方波的频域波形图);grid on;xlabel(频率);ylabel(幅值);t=0:0.001:0.2;N=256;FS=300;w=sawtooth(2*pi*50*t,0.5); figure;plot(t,w); grid on;title(三角波的时域波形图);%产生信号y=fft(w,N); %FFT运算mag=abs(y); %取幅值f=(0:length(y)-1)*FS/length(y);figure;plot(f(1:N/2),mag(1:N/2); %输出FS/2点幅频谱图title(三角波的频域波形图);grid on;xlabel(频率);ylabel(幅值);4) 整周期和非整周期采样时两者的比较，如图19图22所示。图 19图 20图 21图 22此部分的MATLAB代码为：FS=900;SF=10;N=1024;n=0:N-1;t=n/FS;x=80.0*cos(2*3.14*SF*t);figure;plot(t,x);xlabel(t);ylabel(y);title(非整周期时域波形);grid;y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换figure;plot(f(1:N/2),mag(1:N/2);%做频谱图xlabel(频率(Hz);ylabel(幅值);title(非整周期幅频谱图N=1024);grid;FS=1024;SF=10;N=1024;n=0:N-1;t=n/FS;x=80.0*cos(2*3.14*SF*t);figure;plot(t,x);xlabel(t);ylabel(y);title(整周期时域波形);grid;y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换figure;plot(f(1:N/2),mag(1:N/2);%做频谱图xlabel(频率(Hz);ylabel(幅值);title(整周期)幅频谱图N=1024);grid;5) 讨论实验结果通过以上分析，采用整周期采样能够去除频谱中的干扰成分，容易得到准确的频率。二、设计计算倒谱的程序，并对一个模拟信号进行倒谱分析已知模拟信号为： x(n)=120.0*COS(2*3.14*SF*n/FS)（1.0+0.8*COS(2*3.14*SF1*n/FS)）SF1SF,如SF=100，SF1=10。 相关程序如下：clc;clf;Fs=300; %采样频率n=0:1:300;N=128;Sf=20; %信号频率xn=120*cos(2*pi*Sf*n/Fs);%产生波形序列yn=rceps(xn);%倒谱分析subplot(211);plot(n,yn);grid on;title(Matlab自带倒谱函数);xlabel(F