机械故障诊断研讨题1.docx
- 文档编号:8901099
- 上传时间:2023-02-02
- 格式:DOCX
- 页数:21
- 大小:204.46KB
机械故障诊断研讨题1.docx
《机械故障诊断研讨题1.docx》由会员分享,可在线阅读,更多相关《机械故障诊断研讨题1.docx(21页珍藏版)》请在冰豆网上搜索。
机械故障诊断研讨题1
研讨1-信号采集、谱分析、窗函数、整周期采样、谱校正
信号采集过程中一般需要考虑以下几个参数:
信号频率、采样频率、采样长度等,不同
参数的选择对于信号采集的效果会产生直接影响,为了掌握信号采集过程中这些参数对采集
过程及其效果产生的影响,可以通过Matlab或C语言对信号采集与分析处理的过程进行仿
真分析,具体要求如下:
利用Matlab或C语言产生信号:
其中:
频率需要考虑低、中、高,典型的如:
f1=30Hz、f2=400Hz、f3=2000Hz。
每位同
学可按照此思路自己取值,与他人不同;
n(t)为白噪声,均值为零,方差为0.7;
频率、幅值、相位任意设定。
要求每人不相同;
对上述等式进行DFFT处理。
讨论:
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*2100*t+70)+white_noise;(1~4问有小问题,因为上式x中相位应该除以180乘以pi变为角度,第五问纠正过来,而且第五问应该分析三个峰值,我只分析了一个,其余两个,可以谱校正取最大值(即峰值)时分别选取N/2,N/4,N/20从而取到第二第三峰值(即另外两个频率最大值)。
)
1)MATLAB程序结果框图:
通过四个图的比较分析,当采样频率
<2
时,即
,波形无法显示在高频率段的状态,由图一、二可以看出波形出现混叠,而由图四、五可知,当
,即
时,波形在
、
、
三个频率时的幅值高出其他频率,因此对于一个包含多个频率成分的信号进行处理时,要求采样频率
必须高于信号频率成分中最高频率
的两倍,即
这就是采样定理。
2)
采样点数一定的情况下,高频谱比低频的峰值频率带更窄,更可以清晰的看出峰值频率,并且两个峰值频率的混叠情况比低频的要少。
二、采样频率
,采样点数N,采样间隔
:
频率采样间隔
决定了频率分辨力。
越小,分辨力越高。
三、设置
=1,
=1,
=2,n(t)均值为零,方差分别为0.7,1.5,3.2
由图可知,白噪声1的主瓣幅值最高,白噪声3幅值最低,说明白噪声的方差越大,时域信号的相应越准确,噪声对信号的频域分析影响越小。
四、加窗对频谱的影响:
五、整周期采样,谱校正
只有信号的截断长度T为待分析信号周期的整数倍时,才可能使谱线落在f0,获得准确的频谱。
此即为整周期采样,而此信号包含多个频率成分,因此取
、
、
的最小公倍数为T=21000Hz
Matlab输出结果:
MATLAB程序:
第一问:
>figure
(1)%采样频率为200Hz
fs=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+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)%采样频率为1000Hz
fs=1000;
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('采样频率为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(y));
xlabel('频率(Hz)');
ylabel('幅值');
set(gca,'xlim',[0,500]);
figure(3)%采样频率为4300Hz
fs=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');
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,2150]);
figure(4)%采样频率为7000Hz
fs=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*2100*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)%白噪声1
fs=7000;
N=1024;%采样点数
v=0:
N-1;
t=v/fs;%时间序列
white_noise=sqrt(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])
%白噪声2
figure
(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);
subplot(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('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])
>>
第四问:
>>clearall;
figure
(1)%采样频率为7000Hz
fs=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*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);
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=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,'xlim',[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])
>>第五问:
>>clearall
figure
(1)%采样频率为7000Hz
fs=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/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('频率(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=%fn\n',f,A*2,P1/pi*180);
%幅值乘以2,是因为原谱为双边谱,乘以2转化为单边谱
%加汉宁窗
figure
(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('频率(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=%fn\n',f2,A2*2*2,P2/pi*180+90);
%后一个2是对汉宁窗的补偿
[fc,Ac,Pc]=SpecCalibration(A1,A2,A3,f2,P2,deltaf);
fprintf('比例内插法谱校正后的值为nf=%f,A=%f,Phase=%fn\n',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=%fn\n',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,Phase=70.602944n
优化算法做谱校正的值为nf=2099.985920,A=9.011067,Phase=69.852541n
>>子程序
%最优化方法谱校正
function[fc,Ac,Pc]=SpecCalibration_Opt(m,f2,fs,deltaf)
L=length(m);
DFTfun=@(fx)sum(m.*exp(-1j*2*pi*(0:
L-1)*fx/fs))/L;
%设定区间
Al=abs(DFTfun(f2-deltaf));Ar=abs(DFTfun(f2+deltaf));
ifAl>=Ar
fmax=f2;
fmin=f2-deltaf;
else
fmax=f2+deltaf;
fmin=f2;
end
%最优化方法求解
myfun=@(fx)-abs(DFTfun(fx));
fc=fminbnd(myfun,fmin,fmax);
Yf=DFTfun(fc);
Ac=abs(Yf);
Pc=angle(Yf);
End
子程序二:
function[fc,Ac,Pc]=SpecCalibration(A1,A2,A3,f2,P2,deltaf)
ifA3>=A1
r=A3/A2;
%汉宁窗
deltaK=(2*r-1)/(r+1);
else
r=A1/A2;
%汉宁窗
deltaK=(1-2*r)/(r+1);
end
%计算校正值
fc=f2+deltaK*deltaf;
Ac=A2/sinc(deltaK)*(1-deltaK^2);
Pc=P2-pi*deltaK;
end
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 机械 故障诊断 研讨