用MATLAB对信号做频谱分析
1.首先学习下傅里叶变换的东西。学高数的时候老师只是将傅里叶变换简单的说了下，并没有深入的讲解。而现在看来，傅里叶变换似乎是信号处理的方面的重点只是呢，现在就先学习学习傅里叶变换吧。
上面这幅图在知乎一个很著名的关于傅里叶变换的文章中的核心插图，我觉得这幅图很直观的就说明了傅里叶变换的实质。时域上的东西直观的反应到了频域上了，很完美的结合到了一起，233333. 无数正弦波叠加，震荡的叠加的最后结果竟然是方波，同理，任何周期性函数竟然都能拆分为傅里叶级数的形式，这样的简介与优雅，真令人折服。
2.MATLAB对信号做频谱分析
代码：(1)对 f1 = Sa(2t)的频谱分析
1 clear;clc; 2 hold on; 3 R=0.05; 4 t=-1.2:R:1.2; 5 t1 = 2*t; 6 f1=sinc(t1); %Sa函数 7 subplot(1,2,1),plot(t,f1) 8 xlabel(\'t\'),ylabel(\'f1\') 9 axis([-2,2,-0.3,1.2]); %写出Sa函数上下限 10 11 N=1000; 12 k=-N:N; 13 W1=40; 14 W=k*W1/N; 15 F=f1*exp(-j*t\'*W)*R; %f1的傅里叶变换 16 F=real(F); %取F的实部 17 subplot(1,2,2),plot(W,F) 18 xlabel(\'W\'),ylabel(\'F(jw)\')
结果如下图：
(2)对 f2 = u(t+2) - u(t-2)的频谱分析
1 R=0.05; 2 t=-3:R:3; 3 f2=(t>=-2)-(t>=2); 4 subplot(1,2,1),plot(t,f2) 5 grid on; 6 xlabel(\'t\'),ylabel(\'f2\') 7 axis([-3,3,-0.5,1.5]); 8 9 N=1000; k=-N:N; 10 W1=40; 11 W=k*W1/N; 12 F=f2*exp(-j*t\'*W)*R; 13 F=real(F); 14 subplot(1,2,2),plot(W,F) 15 grid on; 16 xlabel(\'W\'),ylabel(\'F(jw)\')
结果如下图：
(3)对f3 = t[u(t+1) - u(t-1) ]的频谱分析
1 R=0.05; 2 h=0.001; 3 t=-1.2:R:1.2; 4 y=t.*(t>=-1)-t.*(t>=1); 5 f4=diff(y)/h; 6 subplot(1,2,1),plot(t,y) 7 xlabel(\'t\'),ylabel(\'y\') 8 axis([-1.2,1.2,-1.2,1.2]); 9 10 N=1000; 11 k=-N:N; 12 W1=40; 13 W=k*W1/N; 14 F=y*exp(-j*t\'*W)*R; 15 F=real(F); 16 subplot(1,2,2),plot(W,F) 17 xlabel(\'W\'),ylabel(\'F(jw)\') 18 axis([-40,40,-0.06,0.06]);
结果如下图：
(4)对正弦波做FFT频谱分析
1 %*************************************************************************% 2 % FFT实践及频谱分析 % 3 %*************************************************************************% 4 %***************正弦波****************% 5 fs=100;%设定采样频率 6 N=128; 7 n=0:N-1; 8 t=n/fs; 9 f0=10;%设定正弦信号频率 10 %生成正弦信号 11 x=sin(2*pi*f0*t); 12 figure(1); 13 subplot(231); 14 plot(t,x);%作正弦信号的时域波形 15 xlabel(\'t\'); 16 ylabel(\'y\'); 17 title(\'正弦信号y=2*pi*10t时域波形\'); 18 grid; 19 20 %进行FFT变换并做频谱图 21 y=fft(x,N);%进行fft变换 22 mag=abs(y);%求幅值 23 f=(0:length(y)-1)\'*fs/length(y);%进行对应的频率转换 24 figure(1); 25 subplot(232); 26 plot(f,mag);%做频谱图 27 axis([0,100,0,80]); 28 xlabel(\'频率(Hz)\'); 29 ylabel(\'幅值\'); 30 title(\'正弦信号y=2*pi*10t幅频谱图N=128\'); 31 grid; 32 33 %求均方根谱 34 sq=abs(y); 35 figure(1); 36 subplot(233); 37 plot(f,sq); 38 xlabel(\'频率(Hz)\'); 39 ylabel(\'均方根谱\'); 40 title(\'正弦信号y=2*pi*10t均方根谱\'); 41 grid; 42 43 %求功率谱 44 power=sq.^2; 45 figure(1); 46 subplot(234); 47 plot(f,power); 48 xlabel(\'频率(Hz)\'); 49 ylabel(\'功率谱\'); 50 title(\'正弦信号y=2*pi*10t功率谱\'); 51 grid; 52 53 %求对数谱 54 ln=log(sq); 55 figure(1); 56 subplot(235); 57 plot(f,ln); 58 xlabel(\'频率(Hz)\'); 59 ylabel(\'对数谱\'); 60 title(\'正弦信号y=2*pi*10t对数谱\'); 61 grid; 62 63 %用IFFT恢复原始信号 64 xifft=ifft(y); 65 magx=real(xifft); 66 ti=[0:length(xifft)-1]/fs; 67 figure(1); 68 subplot(236); 69 plot(ti,magx); 70 xlabel(\'t\'); 71 ylabel(\'y\'); 72 title(\'通过IFFT转换的正弦信号波形\'); 73 grid;
执行结果如下图：
(5)对矩形波做FFT频谱分析
1 %****************2.矩形波****************% 2 fs=10;%设定采样频率 3 t=-5:0.1:5; 4 x=rectpuls(t,2); 5 x=x(1:99); 6 figure(1); 7 subplot(231); plot(t(1:99),x);%作矩形波的时域波形 8 xlabel(\'t\'); 9 ylabel(\'y\'); 10 title(\'矩形波时域波形\'); 11 grid; 12 13 %进行FFT变换并做频谱图 14 y=fft(x);%进行fft变换 15 mag=abs(y);%求幅值 16 f=(0:length(y)-1)\'*fs/length(y);%进行对应的频率转换 17 figure(1); 18 subplot(232); 19 plot(f,mag);%做频谱图 20 xlabel(\'频率(Hz)\'); 21 ylabel(\'幅值\'); 22 title(\'矩形波幅频谱图\'); 23 grid; 24 25 %求均方根谱 26 sq=abs(y); 27 figure(1); 28 subplot(233); 29 plot(f,sq); 30 xlabel(\'频率(Hz)\'); 31 ylabel(\'均方根谱\'); 32 title(\'矩形波均方根谱\'); 33 grid; 34 35 %求功率谱 36 power=sq.^2; 37 figure(1); 38 subplot(234); 39 plot(f,power); 40 xlabel(\'频率(Hz)\'); 41 ylabel(\'功率谱\'); 42 title(\'矩形波功率谱\'); 43 grid; 44 45 %求对数谱 46 ln=log(sq); 47 figure(1); 48 subplot(235); 49 plot(f,ln); 50 xlabel(\'频率(Hz)\'); 51 ylabel(\'对数谱\'); 52 title(\'矩形波对数谱\'); 53 grid; 54 55 %用IFFT恢复原始信号 56 xifft=ifft(y); 57 magx=real(xifft); 58 ti=[0:length(xifft)-1]/fs; 59 figure(1); 60 subplot(236); 61 plot(ti,magx); 62 xlabel(\'t\'); 63 ylabel(\'y\'); 64 title(\'通过IFFT转换的矩形波波形\'); 65 grid;
执行结果如下图：
(6)对白噪声做频谱分析
1 %****************3.白噪声****************% 2 fs=10;%设定采样频率 3 t=-5:0.1:5; 4 x=zeros(1,100); 5 x(50)=100000; 6 figure(1); 7 subplot(231); 8 plot(t(1:100),x);%作白噪声的时域波形 9 xlabel(\'t\'); 10 ylabel(\'y\'); 11 title(\'白噪声时域波形\'); 12 grid; 13 14 %进行FFT变换并做频谱图 15 y=fft(x); %进行fft变换 16 mag=abs(y);%求幅值 17 f=(0:length(y)-1)\'*fs/length(y);%进行对应的频率转换 18 figure(1); 19 subplot(232); 20 plot(f,mag);%做频谱图 21 xlabel(\'频率(Hz)\'); 22 ylabel(\'幅值\'); 23 title(\'白噪声幅频谱图\'); 24 grid; 25 26 %求均方根谱 27 sq=abs(y); 28 figure(1); 29 subplot(233); 30 plot(f,sq); 31 xlabel(\'频率(Hz)\'); 32 ylabel(\'均方根谱\'); 33 title(\'白噪声均方根谱\'); 34 grid; 35 36 %求功率谱 37 power=sq.^2; 38 figure(1); 39 subplot(234); 40 plot(f,power); 41 xlabel(\'频率(Hz)\'); 42 ylabel(\'功率谱\'); 43 title(\'白噪声功率谱\'); 44 grid; 45 46 %求对数谱 47 ln=log(sq); 48 figure(1); 49 subplot(235); 50 plot(f,ln); 51 xlabel(\'频率(Hz)\'); 52 ylabel(\'对数谱\'); 53 title(\'白噪声对数谱\'); 54 grid; 55 56 %用IFFT恢复原始信号 57 xifft=ifft(y); 58 magx=real(xifft); 59 ti=[0:length(xifft)-1]/fs; 60 figure(1); 61 subplot(236); 62 plot(ti,magx); 63 xlabel(\'t\'); 64 ylabel(\'y\'); 65 title(\'通过IFFT转换的白噪声波形\'); 66 grid;
执行结果如下：
请发表评论