%% ======================================================================
%   实验内容 1：方波信号的分解与合成（合集脚本）
%   包括：
%   (1) 方波谐波分解（实验 1.1）
%   (2) 方波谐波合成对比（实验 1.2）
% ======================================================================

clc; clear; close all;

%% ----------------------------------------------------------------------
%   实验内容 1.1：方波信号的傅里叶级数分解
% ----------------------------------------------------------------------

A = 3;                % 幅值
f = 50;               % 基频 50Hz
w0 = 2*pi*f;          % 角频率
T = 1/f;

t = 0 : 0.0001 : 2*T; % 时间：两个周期

harmonics_idx = [1,3,5,7,9];
N_terms = length(harmonics_idx);
harmonics_data = zeros(N_terms, length(t));

coeff_base = 4*A/pi;

% 计算谐波分量
for i = 1:N_terms
    k = harmonics_idx(i);
    harmonics_data(i,:) = (coeff_base / k) * sin(k*w0*t);
end

% ------------------------ 图 1：各谐波分量（分图显示） ------------------
figure('Name','方波分解 - 单独谐波','Color','w');
for i = 1:N_terms
    subplot(N_terms,1,i);
    plot(t, harmonics_data(i,:), 'LineWidth',1.5);
    grid on;
    ylabel(['k=',num2str(harmonics_idx(i))]);
    if i==1
        title('方波的前 5 项谐波分量（50Hz，幅值 3）');
    end
end
xlabel('时间 t (s)');

% ------------------------ 图 2：谐波的三维 mesh 图 ----------------------
figure('Name','方波分解 - 3D Mesh','Color','w');
[X, Y] = meshgrid(t, harmonics_idx);
mesh(X, Y, harmonics_data);
view(3);
title('方波各次谐波的三维波形');
xlabel('时间 t (s)');
ylabel('谐波序号 k');
zlabel('幅值');
colormap jet;

% ------------------------ 图 3：谐波叠加 2D ------------------------------
figure('Name','方波分解 - 2D叠加','Color','w');
hold on;
colors = ['r','g','b','m','k'];
for i=1:N_terms
    plot(t, harmonics_data(i,:), 'Color',colors(i),'LineWidth',1.3);
end
title('前五谐波在同一坐标系叠加');
xlabel('时间 t (s)');
ylabel('幅值');
legend("k=1","k=3","k=5","k=7","k=9");
grid on;


%% ----------------------------------------------------------------------
%   实验内容 1.2：方波信号的傅里叶级数合成
% ----------------------------------------------------------------------
% =========================== 执行两种情况 ===============================

% x1(t): |t|<1 → 1 , |t|>1 → 0
% tau = 2, T = 4
synth_square_wave(4, 2, [5,20,100], 'x1(t): T=4, 宽度 2');

% x2(t): |t|<0.5 → 1 , |t|>0.5 → 0
% tau = 1, T = 2
synth_square_wave(2, 1, [5,20,100], 'x2(t): T=2, 宽度 1');

% =================== 定义一个用于合成的函数 ============================
function synth_square_wave(T, tau, N_values, fig_title)

    t_limit = 1.5 * T;
    t = -t_limit : 0.001 : t_limit;
    w0 = 2*pi/T;

    % 构造理想方波
    original = zeros(size(t));
    for i = 1:length(t)
        t_mod = mod(t(i) + T/2, T) - T/2;
        if abs(t_mod) < tau/2
            original(i) = 1;
        end
    end

    figure('Name', fig_title, 'Color', 'w');

    for k = 1:length(N_values)
        N = N_values(k);
        f_synth = zeros(size(t));

        % 直流项
        a0 = tau/T;
        f_synth = f_synth + a0;

        % 余弦系数 an
        for n = 1:N
            an = 2*(tau/T)*sinc(n*tau/T);
            f_synth = f_synth + an*cos(n*w0*t);
        end

        subplot(length(N_values),1,k);
        plot(t, original,'k--','LineWidth',1); hold on;
        plot(t, f_synth,'r','LineWidth',1.5);
        title(['谐波项数 N = ',num2str(N)]);
        legend('理想方波','合成波形');
        ylim([-0.2,1.2]);
        grid on;
    end
end