matlab信号频谱分析.docVIP

clf; ymax=max(a1); ymin=min(a1); ymean=mean(a1); ystd=std(a1); y=[ymax ymin ymean ystd]; R=1024; a2=[0 0 a1]; a3=[0 a1 0]; a4=[a1 0 0]; a5=(a2+a3+a4)/3; subplot(2,2,1); plot(a0,a1,b-.,a0,a5(2:R+1),r--); xlabel(红色平滑); subplot(2,2,2); m=fft(a1,1024); plot(abs(m)); xlabel(幅值谱); fs=2400; t=0:1/fs:2; subplot(2,2,3); pwelch(a1,60,30,[ ],fs,onesided); xlabel(功率谱); lag=100; [c,lags]=xcorr(a1,lag, unbiased); subplot(2,2,4); plot(lags/fs,c); xlabel(自相关); %FFT变换,获得采样数据基本信息，时域图，频域图 fprintf( 平均加速度 = %7.3f m/ss\n,mean(a1)) %输出本次采样被测量平均值 subplot(2,1,1); plot(a0,a1) %显示原始数据曲线图 axis([min(a0) max(a0) 1.1*floor(min(a1)) 1.1*ceil(max(a1))]) %优化坐标，可有可无 xlabel(时间 (s)); ylabel(加速度（m/ss）); title(振动信号(时域)); grid on; %傅立叶变换 Fs=1024; %得到原始数据时，仪器的采样频率。其实就是length(x)/(max(x)-min(x)); N=1024; %data.txt中的被测量个数，即采样个数。其实就是length(y); z=fft(a1); %频谱分析 f=(0:N-1)*Fs/N; Mag=2*abs(z)/N; %幅值，单位同被测变量y Pyy=Mag.^2; %能量；对实数系列X，有 X.*X=X.*conj(X)=abs(X).^2=X.^2， %显示频谱图(频域) subplot(2,1,2) plot(f(1:N/2), Pyy(1:N/2),r) %显示频谱图 %将这里的Pyy改成Mag就是 幅值-频率图了 axis([min(f(1:N/2)) max(f(1:N/2)) 1.1*floor(min(Pyy(1:N/2))) 1.1*ceil(max(Pyy(1:N/2)))]) xlabel(频率 (Hz)); ylabel(能量); title(频谱图(频域)); grid on; %返回最大能量对应的频率和周期值 [a b]=max(Pyy(1:N/2)); fprintf(\n傅立叶变换结果：\n) fprintf( FFT_f = %1.3f Hz\n,f(b)) %输出最大值对应的频率 fprintf( FFT_T = %1.3f s\n,1/f(b)) %输出最大值对应的周期 1.轴坐标转换 f=(0:length(y)-1)*fs/length(y);%进行对应的频率转换 plot(f,mag);%做频谱图 %****************2.矩形波****************% fs=10;%设定采样频率 t=-5:0.1:5; x=rectpuls(t,2); x=x(1:99); figure(2); subplot(231); plot(t(1:99),x);%作矩形波的时域波形 xlabel(t); ylabel(y); title(矩形波时域波形); grid; %进行FFT变换并做频谱图 y=fft(x);%进行fft变换 mag=abs(y);%求幅值 f=(0:length(y)-1)*fs/length(y);%进行对应的频率转换 figure(2); subplot(232); plot(f,mag);%做频谱图 xlabel(频率(Hz)); ylabel(幅值); title(矩形波幅频谱图); grid; %求均方根谱