机械故障诊断--研讨题1

1、研讨1-信号采集、谱分析、窗函数、整周期采样、谱校正信号采集过程中一般需要考虑以下几个参数：信号频率、采样频率、采样长度等，不同参数的选择对于信号采集的效果会产生直接影响，为了掌握信号采集过程中这些参数对采集过程及其效果产生的影响，可以通过Matlab 或C 语言对信号采集与分析处理的过程进行仿真分析，具体要求如下：利用Matlab 或C 语言产生信号：其中：频率需要考虑低、中、高，典型的如：f1=30Hz、 f2=400Hz、f3=2000Hz。每位同学可按照此思路自己取值，与他人不同；n(t) 为白噪声，均值为零，方差为0.7；频率、幅值、相位任意设定。要求每人不相同；对上述等式进行DFF

2、T 处理。讨论：1）通过设置不同的采样频率，画出时域波形和傅里叶变换后的幅频谱图，用数据分析验证采样定理。讨论在采样点数一定（2 的整数次方）的情况下，如1024 点、2048 点、4096点，采样频率对信号时域复现、频域分析的影响；2）采样频率、采样长度（采样点数）与频率分辨率的关系；3）通过设置不同幅值的信号与噪声，讨论噪声对信号时域分析和频域分析的影响；4）考虑矩形窗和汉宁窗对频谱的影响；5）整周期采样及谱校正。一、 分别选择采样频率200、1000、4300、7000Hz，选择x=4*sin(2*pi*40*t+10)+7*sin(2*pi*500*t+50)+9*sin(2*pi*2

3、100*t+70)+white_noise;(14问有小问题，因为上式x中相位应该除以180乘以pi变为角度，第五问纠正过来，而且第五问应该分析三个峰值，我只分析了一个，其余两个，可以谱校正取最大值（即峰值）时分别选取N/2，N/4,N/20从而取到第二第三峰值（即另外两个频率最大值）。)1) MATLAB程序结果框图:通过四个图的比较分析，当采样频率 figure(1) %采样频率为200Hzfs=200; N=1024;%采样点数v=0:N-1;t=v/fs;%时间序列 white_noise=0+sqrt(0.7)*randn(1,length(t); x=4*sin(2*pi*40*t

4、+10)+7*sin(2*pi*500*t+50)+9*sin(2*pi*2100*t+70)+white_noise; subplot(211);plot(t,x);title(采样频率为200Hz);xlabel(t);ylabel(y);number=1024;y=fft(x,number)/number;n=0:length(y)-1;f=fs*n/length(y);subplot(212);plot(f,abs(y);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,100);figure(2) %采样频率为1000Hzfs=1000; N=1024;

5、%采样点数v=0:N-1;t=v/fs;%时间序列 white_noise=0+sqrt(0.7)*randn(1,length(t); x=4*sin(2*pi*40*t+10)+7*sin(2*pi*500*t+50)+9*sin(2*pi*2100*t+70)+white_noise; subplot(211);plot(t,x);title(采样频率为1000Hz);xlabel(t);ylabel(y);number=1024;y=fft(x,number)/number;n=0:length(y)-1;f=fs*n/length(y);subplot(212);plot(f,abs

6、(y);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,500);figure(3) %采样频率为4300Hzfs=4300; N=1024;%采样点数v=0:N-1;t=v/fs;%时间序列 white_noise=0+sqrt(0.7)*randn(1,length(t); x=4*sin(2*pi*40*t+10)+7*sin(2*pi*500*t+50)+9*sin(2*pi*2100*t+70)+white_noise; subplot(211);plot(t,x);title(采样频率为4300Hz);xlabel(t);ylabel(y);numb

7、er=1024;y=fft(x,number)/number;n=0:length(y)-1;f=fs*n/length(y);subplot(212);plot(f,abs(y);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,2150);figure(4) %采样频率为7000Hzfs=7000; N=1024;%采样点数v=0:N-1;t=v/fs;%时间序列 white_noise=0+sqrt(0.7)*randn(1,length(t); x=4*sin(2*pi*40*t+10)+7*sin(2*pi*500*t+50)+9*sin(2*pi*21

8、00*t+70)+white_noise; subplot(211);plot(t,x);title(采样频率为7000Hz);xlabel(t);ylabel(y);number=1024;y=fft(x,number)/number;n=0:length(y)-1;f=fs*n/length(y);subplot(212);plot(f,abs(y);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,3500);第三问： figure(1) %白噪声1fs=7000; N=1024;%采样点数v=0:N-1;t=v/fs;%时间序列 white_noise=s

9、qrt(0.7)*randn(1,length(t); x=1*sin(2*pi*40*t+10)+1*sin(2*pi*500*t+50)+2*sin(2*pi*2100*t+70)+white_noise; subplot(211);plot(t,x);title(白噪声1);xlabel(t);ylabel(y);number=1024;y=fft(x,number)/number;n=0:length(y)-1;f=fs*n/length(y);subplot(212);plot(f,abs(y);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,fs/2

10、) %白噪声2figure(2)fs=7000;N=1024;%采样点数v=0:N-1;t=v/fs;%时间序列 white_noise=sqrt(1.5)*randn(1,length(t); x=1*sin(2*pi*40*t+10)+1*sin(2*pi*500*t+50)+2*sin(2*pi*2100*t+70)+white_noise; subplot(211);plot(t,x);title(白噪声2);xlabel(t);ylabel(y);number=1024;y=fft(x,number)/number;n=0:length(y)-1;f=fs*n/length(y);s

11、ubplot(212);plot(f,abs(y);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,fs/2) %白噪声3 figure(3)fs=7000;N=1024;%采样点数v=0:N-1;t=v/fs;%时间序列 white_noise=0+sqrt(3.2)*randn(1,length(t); x=1*sin(2*pi*40*t+10)+1*sin(2*pi*500*t+50)+2*sin(2*pi*2100*t+70)+white_noise; subplot(211);plot(t,x);title(白噪声3);xlabel(t);ylabel

12、(y);number=1024;y=fft(x,number)/number;n=0:length(y)-1;f=fs*n/length(y);subplot(212);plot(f,abs(y);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,fs/2)第四问： clear all;figure(1) %采样频率为7000Hzfs=7000;%采样频率N=1024;%采样点数v=0:N-1;t=v/fs;%时间序列 white_noise=sqrt(0.7)*randn(1,length(t); x=4*sin(2*pi*40*t+10)+7*sin(2*pi

13、*500*t+50)+9*sin(2*pi*2100*t+70)+white_noise; subplot(321);plot(t,x);xlabel(t);ylabel(y);title(原始信号);number=1024;y1=(fft(x,number)/number;n=0:length(y1)-1;f1=fs*n/length(y1);subplot(322);plot(f1,abs(y1);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,fs/2)windowB=boxcar(length(x);subplot(323);plot(windowB);

14、xlabel(t);ylabel(y);subplot(324);window1=fft(windowB,number)/number;f2=fs*n/length(window1);plot(f2,abs(window1);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,fs/2)windowB_x=windowB.*x;subplot(325);plot(t,windowB_x);title(加矩形窗后信号);xlabel(t);ylabel(y);subplot(326);fwindowB_x=fft(windowB_x,number)/number;f3

15、=fs*n/length(fwindowB_x);plot(f3,abs(fwindowB_x);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,fs/2)figure(2)windowH=hann(length(x);subplot(221);plot(windowH);xlabel(t);ylabel(y);subplot(222);window2=fft(windowH,number)/number;f4=fs*n/length(window2);plot(f4,abs(window2);xlabel(频率(Hz);ylabel(幅值);set(gca,x

16、lim,0,fs/2)windowH_x=windowH.*x;subplot(223);plot(t,windowH_x);title(加汉宁窗后信号);xlabel(t);ylabel(y);subplot(224);fwindowH_x=fft(windowH_x,number)/number;f5=fs*n/length(fwindowH_x);h1=abs(fwindowH_x);plot(f5,h1);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,fs/2)第五问： clear allfigure(1) %采样频率为7000Hzfs=7000;%采

17、样频率N=1024;%采样点数v=0:N-1;t=v/fs;%时间序列 white_noise=sqrt(0.7)*randn(1,length(t);x=4*sin(2*pi*40*t+10/180*pi)+7*sin(2*pi*500*t+50/180*pi)+9*sin(2*pi*2100*t+70/180*pi)+white_noise; title(原始信号);number=N;YF=(fft(x,number)/number;n=0:length(YF)-1;f1=fs*n/length(YF);subplot(211);Af=abs(YF);plot(f1,Af);xlabel(

18、频率(Hz);ylabel(幅值);set(gca,xlim,0,fs/2)Phf=angle(YF);subplot(212);plot(f1,Phf/pi*180);xlabel(频率(Hz);ylabel(相位);set(gca,xlim,0,fs/2)A, freqind0=max(Af(1:floor(N/2);f = f1(freqind0); P1 = Phf(freqind0);fprintf(直接从频谱读值得到nf = %f, A = %f, Phase = %fnn, f, A*2, P1/pi*180);% 幅值乘以2,是因为原谱为双边谱，乘以2转化为单边谱%加汉宁窗fi

19、gure(2)windowB=hanning(length(x);subplot(211);windowB_x=windowB.*x;title(加矩形窗后信号);fwindowB_x=fft(windowB_x,number)/number;f3=fs*n/length(fwindowB_x);h1=abs(fwindowB_x);%h1为加窗后幅值plot(f3,h1);xlabel(频率(Hz);ylabel(幅值);set(gca,xlim,0,fs/2)subplot(212);phf0=angle(fwindowB_x);plot(f3,phf0/pi*180);xlabel(频率

20、(Hz);ylabel(相位);set(gca,xlim,0,fs/2)%比例内插法谱校正A2, freqind1=max(h1(1:floor(N/2);f2 = f3(freqind1); P2 = phf0(freqind1);A1 = h1(freqind1-1); A3 = h1(freqind1+1);deltaf = diff(f3(1:2);fprintf(加汉宁窗后直接从频谱读值得到nf = %f, A = %f, Phase = %fnn, f2, A2*2*2, P2/pi*180+90); % 后一个2是对汉宁窗的补偿fc, Ac, Pc = SpecCalibrati

21、on(A1, A2, A3, f2, P2, deltaf);fprintf(比例内插法谱校正后的值为nf = %f, A = %f, Phase = %fnn, fc, Ac*2*2, Pc/pi*180+90); % 后一个2是对汉宁窗的补偿% 最优化方法补偿m=windowB_x-mean(windowB_x);fc1, Ac1, Pc1 = SpecCalibration_Opt(m, f2, fs, deltaf);fprintf(优化算法做谱校正的值为nf = %f, A = %f, Phase = %fnn, fc1, Ac1*2*2, Pc1/pi*180+90); % 后一个2是对汉宁窗的补偿% 比例内插法谱校正 直接从频谱读值得到nf = 2098.632813, A = 8.455765, Phase = 15.649589n加汉宁窗后直接从频谱读值得到nf = 2098.632813, A = 8.783132, Phase = 105.437155n比例内插法谱校正后的值为nf = 2099.955726, A = 8.998381, Ph