⽤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)')
View Code
结果如下图:
(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)')
View Code
结果如下图:
(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]);
View Code
结果如下图:
(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);
短时傅里叶变换matlab程序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;
View Code
执⾏结果如下图:
(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;
View Code
执⾏结果如下图:
(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);
版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。
发表评论