使用AnyMathDSP对多频信号进行全面的频谱和相关分析

频率分析是一种基本的信号处理工具：它允许您识别谐波分量，评估其幅度和相位，并探索频谱随时间的演变。对于复杂的信号特性，不仅使用经典的频率响应和频率响应，而且还使用诸如群延迟，功率谱密度，相干性，传递函数以及相关性和时频表示的导出参数。

这个例子分析了一个合成信号--三个正弦信号（50、120、200hz；振幅1.0、0.6、0.3；相位0°、45°、60°）和高斯噪声（方差0.05）的总和。采样频率1000赫兹，持续时间1秒。这样一组模拟了有用成分和干扰的典型混合，这使得在接近真实条件下测试频谱方法的有效性成为可能。

计算在AnyMath环境中使用AnyMathDSP和统计库进行。使用加窗加权(Hanna)、PSD周期图估计、频谱图构造(包括具有布莱克曼窗)和自相关计算。结果通过交互式频谱分析仪可视化为PNG图。

这项工作的目的是对频谱和相关方法与信号参数定量评估相结合的实际演示。计算特征列表:

原始信号是时间表示（p1）
振幅谱**（线性标度）**-带峰值标记（p2）
振幅谱（对数标度）-单位为dB（p3）
相位谱-in度(p4)
PSD-以dB/Hz为单位的频谱功率密度(p6)
累积频谱功率是频谱的累积归一化能量(p7)
群延迟-相位频率导数(p8)
相干性是与50Hz参考信号（p9）的线性耦合程度
传递函数（模块） – |H(f)| (p11)
传递函数(相位)-arg[H(f)]以度为单位(p12)
具有汉纳窗口的频谱图是时频表示（p10，存储为 spectrogram.png)
具有Blackman窗口的频谱图与另一个窗口（p14，保存为 spectrogram_blackman.png)
自相关函数是信号与自身的相关性（p13，存储为 autocorrelation.png)

using EngeeDSP, Statistics

接下来，让我们继续示例的实现，并分析下面介绍的代码块。第一阶段是测试信号的形成，下面的代码设置1000Hz的采样频率，持续时间为1s，形成时间轴。确定三个幅度为1.0、0.6、0.3和相位为0°、45°、60°的正弦曲线（50、120、200Hz），对其进行汇总并添加方差为0.05的高斯噪声。绘制随时间变化的接收信号。

Fs = 1000.0          # 采样率
T = 1.0              # 持续时间
N = Int(Fs * T)      # 计数数目
t = range(0, T - 1/Fs, length=N)  # 时间轴

f1, f2, f3 = 50.0, 120.0, 200.0
A1, A2, A3 = 1.0, 0.6, 0.3
phi1, phi2, phi3 = 0.0, π/4, π/3

signal = A1*sin.(2π*f1*t .+ phi1) +
         A2*sin.(2π*f2*t .+ phi2) +
         A3*sin.(2π*f3*t .+ phi3) +
         0.05*randn(N)

p1 = plot(t, signal,
          title="原始信号（正弦+噪声之和）",
          xlabel="时间，从", ylabel="的振幅",
          legend=false, linewidth=1.5)

振幅谱(线性刻度)

下面的代码将汉纳窗口应用于信号以减少频谱泄漏，执行前向FFT和零频率居中的频谱移位。生成-500至499Hz的频率轴，通过窗口校正计算振幅谱，并用检测到的峰值的标记在线性标度上绘制频率响应。

幅度-频率响应(afc)是信号频谱分量的幅度对频率的依赖性。它使我们能够量化每个谐波分量的能量贡献，识别主导频率，噪声分量的存在和非线性失真。在信号处理任务中，频率响应是滤波、解调、系统状态诊断（振动分析、声学）和信号源识别的基础。线性尺度上的频率响应的构造提供了绝对振幅的可视化表示，这对于测量路径的校准和与参考值的比较是必要的。

window = EngeeDSP.Functions.hann(N)
signal_windowed = signal .* window

Y = EngeeDSP.Functions.fft(signal_windowed)
Y_shifted = EngeeDSP.Functions.fftshift(Y)

freq = collect(range(-Fs/2, Fs/2 - Fs/N, length=N))
freq_shifted = freq

amplitude = abs.(Y_shifted) * 2 / sum(window)
amplitude[1] = amplitude[1] / 2

pos_idx = findall(freq_shifted .> 0)
pos_amp = amplitude[pos_idx]
pos_freq = freq_shifted[pos_idx]

peaks_mask = falses(length(pos_amp))
for i in 2:length(pos_amp)-1
    if pos_amp[i] > pos_amp[i-1] && pos_amp[i] > pos_amp[i+1] && pos_amp[i] > 0.1 * maximum(pos_amp)
        peaks_mask[i] = true
    end
end
peak_freqs = pos_freq[peaks_mask]
peak_amps = pos_amp[peaks_mask]

p2 = plot(freq_shifted, amplitude,
          title="振幅谱",
          xlabel="频率，赫兹", ylabel="的振幅",
          xlims=(0, Fs/2), legend=false,
          linewidth=1.5, color=:blue)
scatter!(p2, peak_freqs, peak_amps, markersize=5, color=:red, label="山峰")

振幅谱（对数标度）和相位谱

下面的代码通过FFT的虚部和实部的反正切计算以弧度为单位的相位谱，并且还通过添加少量ε将幅度谱转换为对数标度（dB），以避免对数从零开始。绘制正频率（从1Hz开始）的频率响应（以分贝为单位）。

频率响应的对数表示（以dB为单位）便于分析具有大动态范围的信号，允许同时观察强和弱谐波分量，以及估计谱密度衰减的速率。 DB量表反映了功率的相对变化（20dB对应于幅度的十倍变化），广泛用于过滤理论，声学和无线电工程来计算传输和衰减系数。

phase = EngeeDSP.Functions.angle.(Y_shifted)
amp_db = 20*log10.(amplitude .+ eps(Float64))
phase_deg = rad2deg.(phase)

p3 = plot(freq_shifted, amp_db,
          title="振幅谱（对数标度）",
          xlabel="频率，赫兹", ylabel="振幅，dB",
          xlims=(1, Fs/2), legend=false,
          linewidth=1.5, color=:green)

相位谱

下面的代码使用先前计算的相位（以弧度为单位）绘制正频率的相位频谱（fch）（以度为单位）。

相位-频率响应(fch)显示频谱分量的初始相位对频率的依赖性。它显示了每个谐波相对于倒计时开始的时间偏移量。频率响应对于评估信号的时间结构至关重要：频率响应的线性度意味着恒定的群延迟和不存在色散失真，这对于通信和信号处理系统非常重要。频率响应图补充了频率响应，给出了复平面中信号的完整频率表示。

p4 = plot(freq_shifted, phase_deg,
          title="相位谱",
          xlabel="频率，赫兹", ylabel="相位，度",
          xlims=(0, Fs/2), legend=false,
          linewidth=1.5, color=:red)

功率谱密度(PSD)

下面的代码使用periodogram方法计算频谱功率密度（PSD）-FFT模的平方，由采样频率和样本数（Fs*N）的乘积归一化。然后，它执行频谱移位以使零频率居中，并将PSD转换为相对于1Hz（dB/Hz）的分贝。构建正频率的PSD图。

频谱功率密度(PSD)表征信号功率在频率上的分布，并确定总功率中有多少落在单个频段(1Hz)上。与幅度谱不同，PSD提供了基于幅度平方的能量解释，这对于具有宽带噪声的随机过程和信号尤其重要。 PSD以dB/Hz为单位的对数表示便于在宽动态范围内可视化。

psd = abs2.(Y) / (Fs * sum(window.^2))
psd_shifted = EngeeDSP.Functions.fftshift(psd)
psd_db = 10*log10.(psd_shifted .+ eps(Float64))

p5 = plot(freq_shifted, psd_db,
          title="光谱功率密度(PSD)",
          xlabel="频率，赫兹", ylabel="PSD, дБ/Гц",
          xlims=(0, Fs/2), legend=false,
          linewidth=1.5, color=:purple)

累积频谱功率

下面的代码将累积频谱功率计算为移位PSD的累积和（cumsum）与其总量的比率，这给出了频谱的归一化累积能量。绘制正频率的累积归一化功率对频率的依赖性。

累积频谱功率（积分能量分布曲线）显示总信号功率的多少集中在从零到给定频率的频率范围内。此特性允许您确定包含给定分数能量（例如，90％功率）的频率范围，这对于评估频谱的有效宽度，设计滤波器和分析单个组件的能量贡献非常重要。

pos_idx = findall(freq_shifted .>= 0)
freq_pos = freq_shifted[pos_idx]
psd_pos = psd_shifted[pos_idx]
cumulative_power_pos = cumsum(psd_pos) / sum(psd_pos)
residual_power = 1 .- cumulative_power_pos

p6 = plot(freq_pos, residual_power,
          title="累积频谱功率（残差）",
          xlabel="频率，赫兹", ylabel="额定功率",
          xlims=(0, Fs/2), legend=false,
          linewidth=2, color=:orange)

群延迟

下面的代码将群延迟计算为相位谱在频率上的负导数（有限差分）。为此，取相邻相位值之间的差值并除以相应频率的差值。群延迟的频率轴定义为相邻频率样本的算术平均值。构造了用于正频率的群延迟的曲线图。

群延迟t(ω)=-dφ(ω)/dw表征每个频率下信号通过系统的时间，并确定窄带信号包络的延迟。对于线性相位谱，群延迟是恒定的，这表明没有色散失真。与常数值的偏差表明频率响应的非线性，导致脉冲信号模糊。

phase_unwrapped = EngeeDSP.Functions.unwrap(phase)
group_delay = -diff(phase_unwrapped) ./ diff(freq_shifted)
freq_gd = (freq_shifted[1:end-1] + freq_shifted[2:end]) / 2

p7 = plot(freq_gd, group_delay,
          title="群延迟",
          xlabel="频率，赫兹", ylabel="延迟，与",
          xlims=(0, Fs/2), legend=false,
          linewidth=1.5, color=:brown)

连贯性

下面的代码生成一个频率为50hz的正弦波形式的参考信号，用相同的加窗加权计算其FFT，然后计算原始信号和参考信号之间的相干性，作为互谱模的平方，由两个信号的频谱功率密度的乘积归一化（加入ε以避免零除）。进行频谱移位，并为从0到1范围内的正频率绘制相干图。

相干性（coherence function）是频域中两个信号之间线性关系的归一化度量，取值从0到1。值1表示完全线性关系（信号在给定频率下完全由参考确定），0表示完全缺乏通信。在这种情况下，选择50Hz正弦波作为参考，因此相干性在50Hz时最大，在其他频率下急剧下降。

reference = sin.(2π*f1*t)
seg_len = 256
overlap = 128
nfft = 512
window_coh = EngeeDSP.Functions.hann(seg_len, "periodic")
step = seg_len - overlap
n_segments = Int(floor((length(signal) - seg_len) / step)) + 1
Pxx = zeros(Float64, nfft)
Pyy = zeros(Float64, nfft)
Pxy = zeros(ComplexF64, nfft)
for i in 1:n_segments
    start = (i-1)*step + 1
    seg_signal = signal[start:start+seg_len-1] .* window_coh
    seg_ref   = reference[start:start+seg_len-1] .* window_coh
    X = EngeeDSP.Functions.fft(seg_signal, nfft)
    Y = EngeeDSP.Functions.fft(seg_ref, nfft)
    Pxx .+= abs2.(X)
    Pyy .+= abs2.(Y)
    Pxy .+= X .* conj.(Y)
end
Pxx ./= n_segments
Pyy ./= n_segments
Pxy ./= n_segments
coherence = abs2.(Pxy) ./ (Pxx .* Pyy .+ eps(Float64))
coherence_shifted = EngeeDSP.Functions.fftshift(coherence)
freq_shifted_coh = EngeeDSP.Functions.fftshift( (0:nfft-1) * Fs / nfft )
pos_idx = freq_shifted_coh .>= 0
freq_pos = freq_shifted_coh[pos_idx]
coh_pos = coherence_shifted[pos_idx]
p8 = plot(freq_pos, coh_pos,
          title="相干性（相对于f1=50Hz）",
          xlabel="频率，赫兹", ylabel="连贯性",
          xlims=(0, Fs/2), legend=false,
          linewidth=1.5, color=:cyan)

结合主要图表和保存

下面的代码将九个图（源信号，频率响应，对数频率响应，频率响应，频率响应，PSD，累积功率，群延迟，相干性）组合成一个3x3网格，将其保存到文件中 "frequency_characteristics.png". 频谱图单独存储（稍后将构建）在 "spectrogram.png".

p_combined = plot(p1, p2, p3, p4, p5, p6, p7, p8,
                  layout=(3, 3), size=(1200, 900),
                  titlefontsize=10, legend=false)
savefig(p_combined, "frequency_characteristics.png")

"/user/start/examples/dsp/spectral_analysis_enseedsp/frequency_characteristics.png"

传递函数（模块）

下面的代码计算传递函数为信号的复频谱与参考信号（50Hz正弦波）频谱的比率，执行频谱移位，选择模块并绘制正频率。

传递函数H(f)=Y(f)/X(f)描述了线性系统在频域中对信号的变换。她的模块 |H(f)| 它显示了每个频率分量在通过系统时的幅度变化的多少倍。在此上下文中，参考信号充当输入，而研究的信号充当输出，这使得能够评估由介质或传输路径引入的失真。

H = Pxy ./ (Pyy .+ eps(Float64))
H_shifted = EngeeDSP.Functions.fftshift(H)
H_mag = abs.(H_shifted)

p11 = plot(freq_pos, H_mag[pos_idx],
           title="传递函数（模块）",
           xlabel="频率，赫兹", ylabel="|H(f)|",
           xlims=(0, Fs/2), legend=false,
           linewidth=1.5, color=:magenta)

传递函数（相位）

下面的代码计算传递函数的相位（频谱复比的角度），将其从弧度转换为度，并绘制正频率的相位响应。

传递函数arg[H(f)]的相位定义了系统（或传播信道）在每个频率下在输入和输出信号之间引入的相移。理想情况下，在没有失真的情况下，相位应该是频率的线性函数，其对应于恒定的群延迟。传递函数的模块和相位的联合分析充分描述了线性系统的频率属性。

H_phase = EngeeDSP.Functions.angle.(H_shifted)

p12 = plot(freq_pos, rad2deg.(H_phase[pos_idx]),
           title="传递函数(相位)",
           xlabel="频率，赫兹", ylabel="阶段H(f),度",
           xlims=(0, Fs/2), legend=false,
           linewidth=1.5, color=:darkgreen)

结合传递函数图

下面的代码将模块的图形和传递函数的相位组合成一个垂直复合图形，并将其保存到一个文件中 "transfer_function.png".

p_combined2 = plot(p11, p12, layout=(2, 1), size=(800, 600))
savefig(p_combined2, "transfer_function.png")

"/user/start/examples/dsp/spectral_analysis_enseedsp/transfer_function.png"

频谱图（汉娜窗口）

下面的代码定义了一个短时傅立叶变换(STFT)函数，该函数具有128个样本的汉纳窗口、64个样本的重叠(50%)和256个FFT长度。该函数将信号分割成重叠的段，应用加窗加权，并计算每个段的FFT，返回频谱矩阵和时间轴。然后对初始信号执行STFT，形成频率轴（0-500Hz），将幅度转换为分贝，并以带有viridis colormap和动态范围限制从-60到10dB的热图的形式构建频谱图。

频谱图（时频表示）显示信号频谱随时间的变化，这对于分析非平稳过程（语音，音乐，振动，瞬态）至关重要。 Hanna窗口通过平滑段边缘来减少频谱泄漏;段重叠（50％）增加时间分辨率并平滑掉伪影;FFT的长度决定了频率分辨率。

window_len = 128                # 窗口长度（定义变量）
noverlap = 64                   # 重叠
nfft = 256                      # FFT长度
window_stft = EngeeDSP.Functions.hann(window_len, "periodic")  # 周期性窗口

function stft(x, window, noverlap, nfft, Fs)
    step = length(window) - noverlap
    n_segments = Int(floor((length(x) - length(window)) / step)) + 1
    S = zeros(ComplexF64, nfft, n_segments)
    t_stft = zeros(n_segments)
    for i in 1:n_segments
        start_idx = (i-1)*step + 1
        end_idx = start_idx + length(window) - 1
        segment = x[start_idx:end_idx] .* window
        S[:, i] = EngeeDSP.Functions.fft(segment, nfft)
        t_stft[i] = (start_idx + length(window)/2) / Fs
    end
    return S, t_stft
end
S, t_stft = stft(signal, window_stft, noverlap, nfft, Fs)
freq_stft = collect(range(0, Fs/2, length=nfft÷2 + 1))
S_db = 20*log10.(abs.(S[1:nfft÷2+1, :]) .+ eps(Float64))

p10 = heatmap(t_stft, freq_stft, S_db,
              title="信号频谱图",
              xlabel="时间，从", ylabel="频率，赫兹",
              color=:viridis, clims=(-60, 10))
savefig(p10, "spectrogram.png")
display(p10)

自相关函数

下面的代码使用xcorr（最大延迟为N–1）计算信号的自相关函数，以秒为单位形成延迟轴，绘制自相关模块并将其保存到文件中 "autocorrelation.png".

自相关函数R(τ)=△s(t)s(t–τ)dt表征时移为τ的信号与自身副本的相似程度。它在零延迟时的最大值等于信号能量;对于周期性信号，自相关保留了周期性，这使得可以在噪声背景下检测隐藏谐波。

acf = EngeeDSP.Functions.xcorr(signal, signal, N-1)
lags = collect(range(-(N-1), N-1, length=2N-1)) / Fs

p13 = plot(lags, abs.(acf),
           title="自相关函数",
           xlabel="延迟，与", ylabel="相关性",
           legend=false, linewidth=1.5, color=:gray)
savefig(p13, "autocorrelation.png")
display(p13)

频谱图（布莱克曼窗口）

下面的代码计算了一个具有256个样本的Blackman窗口，128个样本的重叠（50％）和512个FFT的频谱图，然后在动态范围为-80的inferno colormap中构建热图。..10dB并将其保存到文件中 "spectrogram_blackman.png".

与Hanna窗口相比，Blackman窗口提供了更低的旁瓣电平（-74dB），这改善了频谱泄漏抑制，但扩展了主瓣，降低了频率分辨率。增加窗口长度（256与128）会增加频率分辨率，但会降低时间分辨率。将频谱图与不同的窗口进行比较，可以为特定信号选择频率和时间分辨率之间的最佳折衷。

window_stft2 = EngeeDSP.Functions.blackman(256, "periodic")
S2, t_stft2 = stft(signal, window_stft2, 128, 512, Fs)
freq_stft2 = collect(range(0, Fs/2, length=512÷2 + 1))
S2_db = 20*log10.(abs.(S2[1:512÷2+1, :]) .+ eps(Float64))

p14 = heatmap(t_stft2, freq_stft2, S2_db,
              title="频谱图（布莱克曼窗口，256个样本）",
              xlabel="时间，从", ylabel="频率，赫兹",
              color=:inferno, clims=(-80, 10))
savefig(p14, "spectrogram_blackman.png")
display(p14)

识别峰值并显示统计数据

下面的代码突出显示正频率，找到幅度谱的局部最大值（峰值），阈值为最大值的10％，输出信号参数（采样频率，持续时间，样本数），检测到的频率分量及其幅度，以

通过搜索频谱峰值，您可以自动识别主要谐波分量、其频率和幅度，这是信号识别和诊断所必需的。阈值滤波器（最大值的0.1）可切断噪声排放。 RMS值表征信号的能量，峰值因子（最大值与RMS之比）是其峰值，这对于评估动态范围和放大路径的要求很重要。

println("Частота дискретизации: $Fs Гц")
println("Длительность: $T с")
println("Количество отсчетов: $N")
println("\nОбнаруженные частотные компоненты:")
for i in 1:length(peak_freqs)
    println("  Частота: $(round(peak_freqs[i], digits=1)) Гц, " *
            "Амплитуда: $(round(peak_amps[i], digits=3))")
end
println("Среднее значение: $(round(mean(signal), digits=4))")
println("Среднеквадратичное (RMS): $(round(EngeeDSP.Functions.rms(signal), digits=4))")
println("Пик-фактор (Crest Factor): $(round(maximum(abs.(signal))/EngeeDSP.Functions.rms(signal), digits=3))")

Частота дискретизации: 1000.0 Гц
Длительность: 1.0 с
Количество отсчетов: 1000

Обнаруженные частотные компоненты:
  Частота: 50.0 Гц, Амплитуда: 0.999
  Частота: 120.0 Гц, Амплитуда: 0.597
  Частота: 200.0 Гц, Амплитуда: 0.302
Среднее значение: -0.0018
Среднеквадратичное (RMS): 0.8506
Пик-фактор (Crest Factor): 2.192

AnyMathDSP。光谱分析仪

下面的代码创建具有以下参数的频谱分析仪对象：采样频率1000Hz，单侧频谱，频谱和频谱图显示，dBm单位，滤波器组方法，遗忘系数为0.9的指数平均，Hanna窗，线性..20dB轴，窗口大小1000×600像素。然后，在10次迭代的周期中，它发送一个信号以累积平均，显示分析器窗口并释放对象。

频谱分析仪是一种交互式工具，用于实时可视化信号的频谱功率密度和频谱图。遗忘系数为0.9的指数平均提供了频谱波动的平滑，保留了对信号变化的响应。滤波器组方法提供了与周期图相比更稳定的频谱估计。 dBm单位（相对于1MW的分贝）允许您在1欧姆的参考负载下测量绝对信号功率。

scope = EngeeDSP.spectrumAnalyzer(
    SampleRate = Fs,
    PlotAsTwoSidedSpectrum = false,
    ViewType = "spectrum-and-spectrogram",
    SpectrumUnits = "dBm",
    ReferenceLoad = 1,
    Method = "filter-bank",
    AveragingMethod = "exponential",
    ForgettingFactor = 0.9,
    FrequencyResolutionMethod = "rbw",
    RBWSource = "auto",
    FFTLengthSource = "auto",
    FFTLength = 1024,
    Window = "hann",
    FrequencyScale = "linear",
    YLimits = (-80, 20),
    ColorLimits = (-80, 20),
    TimeResolutionSource = "auto",
    TimeSpanSource = "property",
    TimeSpan = 0.5,
    Size = (1000, 600),
    Title = "Спектр и спектрограмма сигнала",
    ShowGrid = true,
    ShowColorbar = true,
    Colormap = "jet",
    ShowLegend = false,
    ChannelNames = ["Signal"],
    LineWidth = 1.5,
    LineStyle = "-"
)

for i in 1:10
    scope(signal)
end
scope |> display
release!(scope)

结论

频谱相关分析证实存在三个谐波（50,120,200Hz），其振幅与指定谐波匹配，误差小于1％。相位谱的线性度和群延迟的恒定性表明没有色散失真。 PSD和累积功率量化组件的能量贡献，并且在该频率下与50Hz参考的相干性最大。频谱图显示随时间的频率稳定性，自相关证实了信号的周期性结构。所应用的方法有效地恢复了信号参数，适用于工程计算和科学任务。