matlab 求全波整流信号的傅里叶系数

clear all;
clf;%一定要加这行代码，否则可能你之前的代码能运行生成正确图形，改错之后，运行结果还是之前的图形，每次运行前必须先做清理
syms t;
T0 = 1;
N=6;%谐波次数，可以任意指定
m = heaviside(t) - heaviside(t -T0);
x = abs(cos(pi * t)) *m;

subplot(221);
fplot(m*x);

[X1, w1] = fourierseries(x, T0, N);

X = [conj(fliplr(X1(2:N))) X1];
w = [-fliplr(w1(2:N)) w1];

subplot(222);
stem(w, X);

subplot(223);
stem(w, abs(X));

subplot(224);
stem(w, [-angle(X1(2:N))  angle(X1)]);

%计算函数 x 的 N 次谐波的傅里叶级数系数，函数的周期为 T0
%数组 X 存放的是傅里叶系数，也就是幅值
%数组 w 存放的是频率
function [X, w]=fourierseries(x, T0, N)
    syms t;%因为传进来的x函数中包含符号 t，所以函数内部也要定义符号变量 t，否则会报错
    for k = 1: N,
        X(k)=int(x*exp(-j*2*pi*(k-1)*t/T0), t, 0, T0)/T0;
        w(k)=(k-1)*2*pi/T0;
    end
end