matlab求幅频图和相频图--利⽤傅⾥叶级数变换求取会出错问题:绘制信号的频谱和相谱
背景:可以将时域信号转为傅⾥叶级数
Q:求傅⾥叶变化
A:按照定义进⾏求取
Q:如何绘图呢?
教材答案:
Q:能否⽤计算机?A: matlab,代码如下
1clear;
2N=256;%取点数
3fs=N;%采样频率, and sample with entire cycle
4n=0:N-1;t=n/fs;%时间序列
5A=10;% amplitude
6%时域信号处理
7Xt=zeros(1,N);%存储
8T0=1;
9for i=1:N
10    if i<=N/2
11        Xt(i)=A-A/(T0/2)*t(i);
12
13    else
14        Xt(i)= A/(T0/2)*(t(i)-T0/2);
15    end
16end
17
18figure
19subplot(3,1,1);
20plot(t,Xt);
21title('原信号');
22xlabel('time/s');
23ylabel('振幅/V');
24
25%% FFT caculation
26Y=fft(Xt-mean(Xt));%,N+1)/(N+1);%傅⾥叶变换27% remenber to subtract the DC value before FFT. 28Y=Y(1:N/2+1);%由于对称性取频谱前半部分即可29mag=abs(Y);
30%幅值修正
31mag(1)=mag(1)/N;
32mag(2:end-1)=2*mag(2:end-1)/N;
33f=(0:N/2)*fs/N;
34subplot(3,1,2);
35%plot(f,mag,'r');
36 semilogx(log(f),mag,'r');
37xlabel('频率/Hz');
38xlim([0 11])
39ylabel('振幅/V');
40title('频谱图');
41%相频图
42Phase=angle(Y)*180/pi;
43%xiang=rad2deg(xiang);
44subplot(3,1,3);
45plot(f,Phase,'g');
46xlabel('频率/Hz');
matlab求傅里叶变换47ylabel('相位/deg');
48 ylim([-90 90])
49title('相频图');
Q:为什么不能直接⽤求解的傅⾥叶级数来算?A: 我尝试了,结果总不对,主要是相位⾓不对。
1clear;
2N=256*1;%取点数
3fs=256;%采样频率
4n=0:N-1;t=n/fs;%时间序列
5A=10;
6%时域信号处理
7Xt=zeros(1,N);%存储
8sum=0;
9for i=1:N
10    for m=1:2:11  %n取到11 阶信号作近似原信号
11        sum=sum+cos(m*2*pi*t(i))/(m*m);
12        if(m>=11)
13            Xt(i)=sum;%近似信号
14            sum=0;
15        end
16    end
17end
18At=A/2+4*A/(pi*pi)*Xt;
19subplot(3,1,1);
20plot(t,At);
21title('原信号');
22Y=fft(At-mean(At));%傅⾥叶变换
23Y=Y(1:N/2+1);%由于对称性取频谱前半部分即可24mag=abs(Y);
25%幅值修正
26mag(1)=mag(1)/N;
27mag(2:end-1)=2*mag(2:end-1)/N;
28f=(0:N/2)*fs/N;
29subplot(3,1,2);
30plot(f,mag,'r');
31xlabel('频率/Hz');
32xlim([0 10])
33ylabel('振幅');
34title('频谱图');
35%相频图
36xiang=angle(Y);
37xiang=rad2deg(xiang);
38subplot(3,1,3);
39plot(f,xiang,'g');
40xlabel('频率/Hz');
41ylabel('相位');
42title('相频图');

版权声明:本站内容均来自互联网,仅供演示用,请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198,我们将在24小时内删除。