数据采集与信号处理大作业

一、基本内容：基于FFT的功率谱分析程序设计与应用1．基本要求1）对一个人为产生的信号进行采用FFT变换方法进行功率谱分析。已知信号x(n)=120.0*COS(2*3.14*SF*n/FS)式中：n=0,1,2……N-1SF---信号频率FS---采样频率其FFT变换结果X(k)可用下面提供的FFT子程序求出，计算功率谱的公式为：W（k）=2(XR(k)2+XI(k)2)/N式中：k=0,1,2……N/2-1XR(k)---X(k)的实部XI(k)---X(k)的虚部请用VB,VC或C++Builder编译器编程，或采用MATLAB计算，或采用高级语言调用MATLAB计算。处理结果为采用窗口显示时域波形和频域波形。此信号的时域谱，频域谱，功率谱如下图所示：0123456-150-100-50050100150tyx=120.0*cos(2*3.14*SF*t)时域波形数据采集与信号处理作业201020304050607080901000123456x104频率(Hz)幅值x=120.0*cos(2*3.14*SF*t)幅频谱图N=1024010203040506070809010001234567x106频率(Hz)功率谱密度x=120.0*cos(2*3.14*sf*t)功率谱密度数据采集与信号处理作业3其MATLAB代码为：FS=200;SF=10;N=1024;n=0:N-1;t=n/FS;x=120.0*cos(2*3.14*SF*t);subplot(221);plot(t,x);xlabel('t');ylabel('y');title('x=120.0*cos(2*3.14*SF*t)时域波形');grid;y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换subplot(222);plot(f(1:N/2),mag(1:N/2));%做频谱图xlabel('频率(Hz)');ylabel('幅值');title('x=120.0*cos(2*3.14*SF*t)幅频谱图N=1024');grid;Py=2*(y.*conj(y))/N;%计算功率谱密度Pysubplot(223)plot(f(1:N/2),Py(1:N/2));xlabel('频率(Hz)');ylabel('功率谱密度');title('x=120.0*cos(2*3.14*sf*t)功率谱密度');grid;2）对实验所采集的转子振动信号进行频谱分析数据采集与信号处理作业405010015020025030035040045050000.20.40.60.811.21.4频率(Hz)幅值转子振动信号频谱图其MATLAB代码：SF=1000;fid=fopen('E:\张勇论文资料\数据采集与信号处理作业\信号处理大作业\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)信号经过均值化处理或不经过均值化处理的结果比较数据采集与信号处理作业50102030405060708090100010002000300040005000600070008000频率(Hz)功率谱密度含有直流分量的余弦曲线未均值化的功率谱波形图010203040506070809010001000200030004000500060007000频率(Hz)功率谱密度含有直流分量的余弦曲线均值化的功率谱波形图通过以上两个图的分析，我们可以看出均值化处理后的频谱的低频段消失，这就去去除了常规的干扰频谱，如环境噪声等，对我们进行频谱分析有很大作用。数据采集与信号处理作业6其MATLAB代码为Fs=200;%采样频率n=0:1:200;N=1024;SF=10;%信号频率xn=120.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）采用不同窗函数时的谱结果（矩形窗函数,汉宁窗函数，汉明窗）0123456-150-100-50050100150ty汉宁窗时域波形数据采集与信号处理作业70102030405060708090100050010001500频率(Hz)幅值汉宁窗频域特性0123456-150-100-50050100150ty汉明窗时域波形数据采集与信号处理作业801020304050607080901000100020003000400050006000频率(Hz)幅值汉明窗频域特性0123456-150-100-50050100150ty矩形窗时域波形数据采集与信号处理作业901020304050607080901000123456x104频率(Hz)幅值矩形窗频域特性其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('汉宁窗频域特性');数据采集与信号处理作业10grid;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）典型函数的频谱（矩形窗函数,汉宁窗函数，直线，阶跃函数，δ函数，方波，三角波等）数据采集与信号处理作业1105010015020025000.511.52矩形窗函数的时域波形图0501001500100200300矩形窗函数频域波形图频率幅值05010015020025030000.20.40.60.81δ函数的时域波形图05010015000.511.52δ函数的频域波形图频率幅值数据采集与信号处理作业1200.020.040.060.080.10.120.140.160.180.200.51方波的时域波形图050100150050100150方波的频域波形图频率幅值05010015020025030000.20.40.60.81汉宁窗函数的时域波形图050100150050100150汉宁窗函数频域波形图频率幅值数据采集与信号处理作业1305010015020025030000.511.52阶跃函数的时域波形图0501001500100200300阶跃函数的频域波形图频率幅值00.020.040.060.080.10.120.140.160.180.2-1-0.500.51三角波的时域波形图050100150020406080三角波的频域波形图频率幅值此部分MATLAB代码如下：t=0:0.001:0.2;N=256;数据采集与信号处理作业14FS=300;w=boxcar(N);%产生信号subplot(211);plot(w);title('矩形窗函数的时域波形图');axis([0,260,0,2]);gridon;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('频率');ylabel('幅值');t=0:0.001:0.2;N=256;FS=300;w=hanning(N);%产生信号subplot(211);plot(w);title('汉宁窗函数的时域波形图');gridon;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('频率');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);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;数据采集与信号处理作业15N=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('频率');ylabel('幅值');t=0:0.001:0.2;N=256;FS=300;w=zeros(1,N);w(1)=1;%产生信号subplot(211);plot(w);grid;title('δ函数的时域波形图');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('频率');ylabel('幅值');t=0:0.001:0.2;N=256;FS=300;w=square(2*pi*50*t);%产生信号subplot(211);plot(t,w);title('方波的时域波形图');axis([0,0.2,-0.2,1.2]