matlab 求幅频图和相频图--利用傅里叶级数变换求取会出错

问题：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函数，得出错误的结果。原因未知！！