问题:matlab绘制信号的频谱和相谱
背景:傅里叶可以将时域信号转为傅里叶级数
Q: 求傅里叶变化
A: 按照定义进行求取
Q: 如何绘图呢?
教材答案:
Q:能否用计算机?
A: matlab, 代码如下
clear;
N=256;%取点数
fs=N;%采样频率, and sample with entire cycle
n=0:N-1;t=n/fs;%时间序列
A=10;% amplitude
%时域信号处理
Xt=zeros(1,N);%存储
T0=1;
for i=1:N
if i<=N/2
Xt(i)=A-A/(T0/2)*t(i);
else
Xt(i)= A/(T0/2)*(t(i)-T0/2);
end
end
figure
subplot(3,1,1);
plot(t,Xt);
title('原信号');
xlabel('time/s');
ylabel('振幅/V');
%% FFT caculation
Y=fft(Xt-mean(Xt));%,N+1)/(N+1);%傅里叶变换
% remenber to subtract the DC value before FFT.
Y=Y(1:N/2+1);%由于对称性取频谱前半部分即可
mag=abs(Y);
%幅值修正
mag(1)=mag(1)/N;
mag(2:end-1)=2*mag(2:end-1)/N;
f=(0:N/2)*fs/N;
subplot(3,1,2);
%plot(f,mag,'r');
semilogx(log(f),mag,'r');
xlabel('频率/Hz');
xlim([0 11])
ylabel('振幅/V');
title('频谱图');
%相频图
Phase=angle(Y)*180/pi;
%xiang=rad2deg(xiang);
subplot(3,1,3);
plot(f,Phase,'g');
xlabel('频率/Hz');
ylabel('相位/deg');
ylim([-90 90])
title('相频图');
Q: 为什么不能直接用求解的傅里叶级数来算?
A: 我尝试了,结果总不对,主要是相位角不对。
clear;
N=256*1;%取点数
fs=256;%采样频率
n=0:N-1;t=n/fs;%时间序列
A=10;
%时域信号处理
Xt=zeros(1,N);%存储
sum=0;
for i=1:N
for m=1:2:11 %n取到11 阶信号作近似原信号
sum=sum+cos(m*2*pi*t(i))/(m*m);
if(m>=11)
Xt(i)=sum;%近似信号
sum=0;
end
end
end
At=A/2+4*A/(pi*pi)*Xt;
subplot(3,1,1);
plot(t,At);
title('原信号');
Y=fft(At-mean(At));%傅里叶变换
Y=Y(1:N/2+1);%由于对称性取频谱前半部分即可
mag=abs(Y);
%幅值修正
mag(1)=mag(1)/N;
mag(2:end-1)=2*mag(2:end-1)/N;
f=(0:N/2)*fs/N;
subplot(3,1,2);
plot(f,mag,'r');
xlabel('频率/Hz');
xlim([0 10])
ylabel('振幅');
title('频谱图');
%相频图
xiang=angle(Y);
xiang=rad2deg(xiang);
subplot(3,1,3);
plot(f,xiang,'g');
xlabel('频率/Hz');
ylabel('相位');
title('相频图');
相位角错了?? 没找到原因!!
修改相位角求法:
xiang=atan(imag(Y)./real(Y));
xiang(1)=0;
还是不对
参考 Matlab求解周期函数的傅里叶级数以及作频谱图与相位图_xbb0720的博客-CSDN博客_matlab傅里叶级数
修改如下:
傅里叶级数函数---按照定义来求,是最正确的方法
function [ f ] = Tfourierseries( a0, an, bn, m, t )
%T FOURIERSERIES
%
% input the values of a0, an ,bn from fourier seriers experssion
% t time
%m is the number of harmonic frequency
f=a0;
for n=1:m
f=f+an(n)*cos(n*pi.*t);% wt
f=f+bn(n)*sin(n*pi.*t);
end
通过查阅CSDN:更好的 傅里叶级数函数定义参看:一个MATLAB中傅里叶级数变换的函数_张雅琼的博客-CSDN博客_matlab 计算傅里叶级数
程序:
获得时域信号的,用多个谐波信号来合成时域信号。
%求傅里叶变换
N=1000;%取点数
fs=200;%采样频率, and sample with entire cycle
num=0:N-1;t=num/fs;%时间序列
A=10;% amplitude
%tripuls为三角波函数,进行偏移叠加处理可以得到一个类似三角周期函数的图
a0= A/2;
Nf=11;% 谐波的阶次
an=zeros(1,Nf);
bn=zeros(1,Nf)
for i=1 :Nf
if mod(i,2) ==0
an(i)=0;
else
an(i)=1/i^2*4*A/(pi*pi);
end
end
%3次谐波叠加
f3=Tfourierseries(a0, an, bn, 3, t);
%9次谐波叠加
f9=Tfourierseries(a0, an, bn, 9, t);
% %21次谐波叠加
% f21=trifourierseries(a0, an, bn, 21, t);
% %45次谐波叠加
% f45=trifourierseries(a0, an, bn, 45, t);
%级数展开图
subplot(2,3,1);plot(t,f3,'b');grid on
%axis([-6.1,6.1,-0.1,1.1]);
title('展开3项');
xlabel('t');ylabel('f(t)');
subplot(2,3,4);plot(t,f9,'b');grid on
title('展开9项');
% xlabel('t');ylabel('f(t)');
% subplot(2,3,2);plot(t,f,'r',t,f21,'b');grid on
% axis([-6.1,6.1,-0.1,1.1]);title('展开21项');
% xlabel('t');ylabel('f(t)');
% subplot(2,3,5);plot(t,f,'r',t,f45,'b');grid on
% axis([-6,6,-0.1,1.1]);title('展开45项');
% xlabel('t');ylabel('f(t)');
绘制频谱图---按照定义来;
%幅度谱--相位谱
num=0:1:10;
%注意A0需要自己赋值
An=sqrt(an.^2+bn.^2);An(1)=a0;
%phi0同理
phi=atan(-bn./an);phi(1)=0;
subplot(2,3,3);stem(num,An,'b');
grid on; xlim([0 10]);
title('幅度谱');xlabel('n');ylabel('An');
ylim([0 1])
subplot(2,3,6);plot(num,phi,'b');
grid on; axis([-0.1,10.1,-0.1,1.1]);
title('相位谱');xlabel('n');ylabel('ψn');
按照定义可以求出相位,如果用phase函数,得出错误的结果。原因未知!!