信号与系统¶

如你所见，信号与系统一共两次大作业，一个5000字的综述和一个MATLAB编程。着重说后一次

~~第二次大作业多少有些 fly bitch的味道~~，这里会提到一些必要的技术，剩下的则需要读者自行学习。

电力系统谐波分析¶

第二次是电力系统谐波分析，在五一左右发布，期末考完交，中间有两次精确度排名。谐波分析当你们学完前半学期就有点思路了，基本上是 FFT 处理原始信号，但是这时候会有精确度的问题(可以看看书的第 10 和 11 章)，就得加窗，之后还得考虑一些优化办法，例如双峰谱线插值&三峰谱线插值等等。小波变换最好别用，不太可行。

如前所述，这次大作业的核心方法是 fft, 所有的修正方法都会围绕它展开。

实现思路¶

学习路线¶

学期中，推荐大家先去简单学一学MATLAB的基本操作。MATLAB语言和C差别不是很大，注意数组、函数、绘图这几大块的学习。

本次大作业主要用到的函数是:fft, hht, plot

实现思路¶

难点在于 fft 的精确计算，由于频谱泄露等一系列原因，对于信号的精确分析往往要借助一些修正方法，例如三谱线插值、频谱细化等等。

谐波计算¶

谐波计算有如下的办法

三谱线插值¶

具体原理知网有论文讲解，自己查阅。这个方法实现起来并不困难，不同的窗函数有不同的误差，hamming窗求解频率精度较高，blackharrisman窗求解谐波精度较高，可以自己尝试。

频谱细化¶

频谱细化的主流方法包括 czt 变换和 ZoomFFT, 这里可以使用 czt 变换。贴出一份往届代码供参考:

function [f0,A,phase] = FundCalculation(x,fs,window)
    N = length(x);
    % 设置不同的窗函数
    w = hamming(N,'periodic');
    m = window; 
    L = 0:1/m:0.5;
    % 连续z变换频谱细化
    W = czt(w, m + 1, exp(-1i * 2 * pi / (m * N)), 1);
    index = 1:length(L);
    Alpha = abs(W(end - index + 1)) ./ abs(W(index));
    % 加窗并获得基波指标
    X = fft(x.*w);
    Xa = abs(X);
    m0 = find(Xa == max(Xa), 1);
    if Xa(m0 + 1) > Xa(m0 - 1)
        m1 = m0 + 1;
    else
        m1 = m0 - 1;
    end

    % 三次样条插值估计频率误差
    alpha = Xa(m1) / Xa(m0);
    M = spline(Alpha, L, alpha);
    if m1 < m0
        M = -M;
    end
    f1 = (m0 + M - 1) / N;
    f0=f1 * fs;
    %disp(["fundamental frequency:",f0]);
    %幅值和相位估计
    k = 0:N-1;
    W1 = sum(w'.*exp(1i * 2 * pi * M * k / N));
    A = abs(2 * X(m0) / W1 / sqrt(2));
    phase = angle(X(m0)) - angle(W1);
end

噪声查找¶

去噪方法: 经验模态分解、小波去噪、平滑... 比如:

小波去噪

% 小波去噪
denoised_signal = wdenoise(Main0,9,Wavelet='sym5',DenoisingMethod="SURE");
B=smoothdata(denoised_signal,'lowess');

经验模态分解

imf = emd(Main0(:,2));
Main=zeros(2,10000)';
Main(:,2) = Main0(:,2) -imf(:,1);

关于噪声分布，可以考虑使用频率直方图绘制。