信号作业1

数据采集与信号处理1、基于FFT的功率谱分析程序设计与应用1).已知信号x(n)=100.0*COS(2*3.14*SF*n/FS)+50.0*COS(4*3.14*SF*n/FS)，n=0,1,2……N-1式中：N-1.已知信号x(n)=100.0*COS(2*3.14*SF*n/FS)+50.0*COS(4*3.14*SF*n/FS)，n=0,1,2……N-1式中：N---采样点数，必须为2xSF---信号频率，SF=10HzFS---采样频率其FFT变换结果X(k)可用下面提供的FFT子程序求出，计算功率谱的公式为：G(k)=2(XR(k)2+XI(k)2)/N,k=0,1,2……N/2-1式中：XR(k)---X(k)的实部XI(k)---X(k)的虚部采用MATLAB方法来算出信号的时域波形和功率谱密度波形。程序如下：FS=300;SF=10;N=1024;n=0:N-1;t=n/FS;x=100.0*cos(2*3.14*SF*t)+50.0*cos(4*3.14*SF*t);subplot(211);plot(t,x);xlabel('t');ylabel('y');title('x=100.0*cos(2*3.14*SF*t)+50.0*cos(4*3.14*SF*t)时域波形');grid;y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换subplot(212);plot(f(1:N/2),mag(1:N/2));%做频谱图xlabel('频率(Hz)');ylabel('幅值');title('x=100.0*cos(2*3.14*SF*t)+50.0*cos(4*3.14*SF*t)幅频谱图N=1024');grid;y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*FS/length(y);%进行对应的频率转换subplot(212);plot(f(1:N/2),mag(1:N/2));%做频谱图xlabel('频率(Hz)');ylabel('幅值');title('x=100.0*cos(2*3.14*SF*t)+50.0*cos(4*3.14*SF*t)幅频谱图N=1024');grid;Py=2*(y.*conj(y))/N;%计算功率谱密度Pysubplot(211)plot(f(1:N/2),Py(1:N/2));xlabel('频率(Hz)');ylabel('功率谱密度');title('x=100.0*cos(2*3.14*sf*t)+50.0*cos(4*3.14*SF*t)功率谱密度');grid;此信号的时域谱和频域谱如下图所示：2.对一个用A/D数据采集板采集的信号进行频谱分析（1）fanbo_45HZ_1024_1000HZ，45Hz方波发生器输出信号采集结果，采样频率1000Hz，采样点数1024MATLAB程序如下：SF=1000;fid=fopen('D:\数据文件\fanbo_45HZ_1024_1000HZ');[a,N]=fscanf(fid,'%f');fclose(fid);y=fft(a,N);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));xlabel('频率(Hz)');ylabel('幅值');title('方波频谱图');grid;05010015020025030035040045050000.20.40.60.811.21.41.61.82频率(Hz)幅值方波频谱图（2）正弦波发生器输出信号采集结果，采样频率1000Hz，采样点数1024.MATLAB程序如下：SF=1000;fid=fopen('D:\数据文件\sin_45HZ_1024_1000HZ');[a,N]=fscanf(fid,'%f');fclose(fid);y=fft(a,N);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));xlabel('频率(Hz)');ylabel('幅值');title('正弦波频谱图');grid;05010015020025030035040045050000.20.40.60.811.21.4频率(Hz)幅值正弦波频谱图（3）A_1_1024_1000HZ，汽轮机轴向振动位移信号采集结果，采样频率1000Hz，采样点数1024MATLAB程序如下：SF=1000;fid=fopen('D:\数据文件\A_1_1024_1000HZ');[a,N]=fscanf(fid,'%f');fclose(fid);y=fft(a,N);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));xlabel('频率(Hz)');ylabel('幅值');title('频谱图');grid;0501001502002503003504004505000100200300400500600频率(Hz)幅值频谱图由上面的三幅图我们可以看出，谐波由少到多的是：正弦波谐波最少，只有一个；三角波居中间；方波的谐波最为丰富，所以方波发生器又叫“多谐振荡器”。3).讨论：①信号经过均值化处理或不经过均值化处理的结果比较。由于对信号测量的时候很容易引入不必要的直流分量等常规的干扰信号，在做变换的时候也可能产生冲击。从下面两幅图可以看出，经过均值化处理后这种类似的干扰信号被去除掉，使信号处理起来更加简单，这对信号处理具有重要的意义。②采用不同窗函数时的谱结果分析（矩形窗函数,汉宁窗函数，汉明窗）。其MATLAB代码为：FS=500;SF=10;N=516;n=0:N-1;t=n/FS;x=80*cos(2*pi*SF*t);wd=boxcar(N);x=x.*wd';y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*FS/length(y);subplot(211);plot(f(1:N/2),mag(1:N/2));xlabel('频率(Hz)');ylabel('幅值');title('矩形窗的频谱图');grid;y=fft(rectangle,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('矩形窗函数的频域波形');xlabel('频率(Hz)');ylabel('幅值');grid;05010015020025000.511.52x104频率(Hz)幅值矩形窗的频谱图050100150200250177178179180矩形窗函数的频域波形频率(Hz)幅值其MATLAB代码为：t=0:0.001:0.4;N=1024;FS=1000;h=hanning(N);%产生信号subplot(211);plot(h);title('汉宁窗函数的时域波形');xlabel('t');ylabel('幅值');grid;其MATLAB代码为：FS=500;SF=20;N=256;n=0:N-1;t=n/FS;x=80*cos(2*pi*SF*t);wd=hamming(N);x=x.*wd';y=fft(x,N);mag=abs(y);f=(0:length(y)-1)*FS/length(y);subplot(313);plot(f(1:N/2),mag(1:N/2));xlabel('频率(Hz)');ylabel('幅值');title('汉明窗的频谱图');grid;③人为产生典型信号并进行谱分析，自由讨论计算结果（矩形窗函数,汉宁窗函数，直线，δ函数，方波等）。矩形窗函数MATLAB程序如下：t=0:0.001:0.4;N=512;FS=200;rectangle=boxcar(N);%产生信号subplot(221);plot(rectangle);title('矩形窗时域波形');axis([0,260,0,2]);xlabel('t');ylabel('幅值');grid;y=fft(rectangle,N);%FFT运算mag=abs(y);%取幅值f=(0:length(y)-1)*FS/length(y);subplot(222);plot(f(1:N/2),mag(1:N/2));%输出FS/2点幅频谱图title('矩形窗频域波形');xlabel('频率(Hz)');ylabel('幅值');grid;汉宁窗函数MATLAB程序如下：t=0:0.001:0.4;N=1024;FS=1000;h=hanning(N);%产生信号subplot(211);plot(h);title('汉宁窗时域波形');xlabel('t');ylabel('幅值');grid;y=fft(h,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('汉宁窗频域波形');xlabel('频率(Hz)');ylabel('幅值');grid;直线MATLAB程序如下：t=0:0.001:0.4;N=256;FS=500;z=t;%产生信号subplot(231);plot(z);title('直线的时域波形');xlabel('t');ylabel('幅值');grid;y=fft(z,N);%FFT运算mag=abs(y);%取幅值f=(0:length(y)-1)*FS/length(y);subplot(232);plot(f(1:N/2),mag(1:N/2));%输出FS/2点幅频谱图title('直线频域波形图');xlabel('频率(Hz)');ylabel('幅值');grid;函数MATLAB程序如下：t=0:0.001:0.4;N=256;FS=550;w=zeros(1,N);w(1)=1;%产生信号subplot(241);plot(w);title('δ函数的时域波形');xlabel('t');ylabel('幅值');grid;y=fft(w,N);%FFT运算mag=abs(y);%取幅值f=(0:length(y)-1)*FS/length(y);subplot(242);plot(f(1:N/2),mag(1:N/2));%输出FS/2点幅频谱图title('δ函数的频域波形');xlabel('频率(Hz)');ylabel('幅值');grid;方波函数MATLAB程序如下：t=0:0.001:0.4;N=256;FS=200;w=square(2*pi*45*t)%产生信号subplot(251);plot(t,w);title('方波的时域波形');axis([0,0.2,-0.2,1.2]);xlabel('t');ylabel('幅值');grid;y=fft(w,N);%FFT运算mag=abs(y);%取幅值f=(0:length(y)-1)*FS/length(y);subplot(252);plot(f(1:N/2),mag(1:N/2));%输出FS/2点幅频谱图title('方波的频域波形');xlabel('频率(Hz)');ylabel