使用EngeeDSP进行FFT和频谱分析

我们继续使用频谱分析示例和快速傅里叶变换（FFT）函数的应用来熟悉EngeeDSP库的功能。

该示例基于来自先前发布的项目的数据（[FFT]（https://engee.com/community/ru/catalogs/projects/osnovy-tsos-bystroe-preobrazovanie-fure-bpf ），[光谱分析]（https://engee.com/community/ru/catalogs/projects/osnovy-tsos-spektralnyi-analiz ））这是使用Julia插件库执行的，例如 DSP.jl, FFTW.jl, Statistics.jl 等等。在这种情况下，不需要连接这些库，因为它们的功能由EngeeDSP复盖。

下面将讨论该功能的应用特性 **fft**评估测试信号的幅度和相位谱，以及频谱分析的实际任务，以确定录制的音频信号的特性。

FFT函数

关于信号的频率和相位组成的信息可以通过使用傅立叶变换方法从分析的时域切换到频域来获得。作为一个基本的测试信号，我们将产生四个具有不同振幅、基频和初始相位的正弦波的总和。:

fs = 1000;                              # частота дискретизации сигнала
dt = 1/fs;                              # период дискретизации
stoptime = 1;                           # время окончания сигнала
t = 0:dt:stoptime-dt;                   # вектор отсчётов времени
amplitude = [0.58 1.1 0.74 0.25];       # вектор четырёх амплитуд
f0 = [44 130 251 305];                  # вектор частот в Гц
phase_deg = [45 -135 90 60];            # вектор начальных фаз в градусах
phase = deg2rad.(phase_deg);            # вектор начальных фаз в радианах

four_sines = amplitude .* cos.(2*pi*t.*f0 .+ phase);
sum_sines = vec(sum(four_sines, dims = 2));
plot(t, sum_sines, title = "Сумма четырёх синусоид", xlim = (0, 0.15), legend = false,
                    xguide = "Время (с)", yguide = "Амплитуда")

让我们将函数应用于信号 **fft**直接和可视化的结果:

A = EngeeDSP.Functions.fft(sum_sines);
plot(A)

该图是在复平面上绘制的，因为FFT输出端的复矢量包含有关分量幅度和相位的信息。还值得注意的是，在算法的输出处的复数样本的数量对应于信号在其输入处的时域中的样本的数量。

让我们在FFT的输出端显示复数向量的模块:

plot(abs.(A), xguide = "Номер отсчёт БПФ", yguide = "Модуль БПФ", legend = false)

功能 fftshift 从EngeeDSP组成中，它允许您以频谱的通常形式表示FFT输出，在我们的情况下，相对于中心零频率对称。但是，我们还必须将离散频率的矢量沿横坐标轴在从以前 :

B = EngeeDSP.Functions.fftshift(A);
nsamp = length(t);
df = fs/nsamp;
println("Разрешение по частоте = ", df, "  Гц")
freq_vec = -fs/2:df:fs/2 - df;
plot(freq_vec, abs.(B), xguide = "Частота (Гц)", yguide = "Модуль БПФ", legend = false)

Разрешение по частоте = 1.0  Гц

最后，为了获得信号的幅度谱，我们需要考虑FFT函数输出处的样本数量，并且还考虑到真实信号的能量也延伸到负频率的区域，在我们的情况下，这是无信让我们使用线性图表显示正频率范围内信号的振幅谱，确保在合成频率下正确计算正弦曲线的振幅。:

S = 2*B/nsamp;
plot(freq_vec, abs.(S), line=:stem, marker=:circle, xlim=(0, 500),
        xguide = "Частота (Гц)", yguide = "Амплитуда",
        title = "Амплитудный спектр тестового сигнала")

让我们通过函数确定FFT样本的"角度"来构建测试信号的相位谱 angle. 我们将用标记标记我们感兴趣的频率分量，确保相位计算正确。:

idx = f0 .+ Int(fs/2 + 1);
plot(freq_vec, rad2deg.(angle.(S)), line=:stem, xlim=(0, 500),
        xguide = "Частота (Гц)", yguide = "Фаза (град)",
        title = "Фазовый спектр тестового сигнала")
scatter!(freq_vec[idx], rad2deg.(angle.(S[idx])), legend = false)

设置分析任务

让我们把吉他弦的录音作为分析的信号. 主要任务是通过分析信号的频谱来确定正在播放哪个音符。让我们用函数加载信号的数值样本 wavread，这也允许您选择信号的采样频率 fs -在我们的例子中，它等于8000赫兹:

using WAV
audio_data, fs = wavread("guitar.wav");

让我们在时域中计算信号的样本数，计算时间步长，并形成信号的样本的时间向量和数值向量。 sig:

nsamples = length(audio_data);
dt = 1/fs;
t = 0:dt:(nsamples - 1) .* dt;
sig = vec(audio_data);

我们将使用一个额外的功能来收听音频。 audioplayer:

using Base64
function audioplayer(s, fs);
    buf = IOBuffer();
    wavwrite(s, buf; Fs=fs);
    data = base64encode(unsafe_string(pointer(buf.data), buf.size));
    markup = """<audio controls="controls" {autoplay}>
                <source src="data:audio/wav;base64,$data" type="audio/wav" />
                Your browser does not support the audio element.
                </audio>"""
    display("text/html", markup);
end

audioplayer (generic function with 1 method)

让我们听一个振荡的字符串的声音:

audioplayer(sig,fs)

让我们显示时域的波动:

plot(t,sig, title = "Колебания струны во времени", 
        xguide = "Время (с)", yguide = "Амплитуда", legend = false)

在时域上通过周期信号的形状找出基波振荡频率是可能的，但仍然非常困难。

信号频谱的估计

让我们使用上面应用的操作序列来构造音频信号功率谱。我们还将使用FFT方法评估频谱，但这次我们将限制函数输出处的点数。 fft 8000计数。

进一步的步骤重复标准算法-计算复杂的FFT，移位，按级别归一化。但要沿纵坐标轴显示，振幅的平方，即功率，已经被取并转换为分贝。:

nFFT = Int(fs);                         # кол-во отсчётов на выходе БПФ
A = EngeeDSP.Functions.fft(sig, nFFT);
B = EngeeDSP.Functions.fftshift(A);
S = (2*B)/nFFT;
df = fs/nFFT;
freq_vec = -fs/2:df:fs/2-df;
plot(freq_vec, 10log10.(abs.(S.^2)),
        xguide = "Частота (Гц)", yguide = "Мощность (дБ)",
        title = "Спектр мощности аудиосигнала", legend = false)

我们观察到音频信号功率谱范围从-4,000Hz到4,000Hz。

谐波结构在光谱上清晰可见-等距峰。谐波是正弦振荡，其频率是基频的倍数。 .

-在信号的最低频率（一次谐波），和高次谐波是频率的形式，在哪里（第二，第三等。).

基频的确定

通过分析谐波，可以确定信号频谱中的基频。但首先，让我们根据使用的算法编写一个自定义频谱估计函数。作为额外的输入参数，它将接受信号的采样频率，FFT样本数和"单向"频谱标志。:

function fftspectrum(sig, fs; nfft = nothing, onesided = nothing)
    if nfft != nothing
        nsamples = nfft;
        sig = sig[1:nfft];
    else
        nsamples = length(sig);
    end
    A = EngeeDSP.Functions.fft(sig);
    B = EngeeDSP.Functions.fftshift(A);
    S = (2*B)/nsamples;
    out = 10log10.(abs.(S.^2));
    df = fs/nsamples;
        if onesided == true
            f = 0:df:fs/2;
            out = out[Int(ceil(nsamples/2)):end]
        else            
            f = -fs/2:df:fs/2-df;
        end
    return out, f
end

fftspectrum (generic function with 1 method)

让我们分析低频范围（从0到1500Hz）。我们将在频谱中找到峰值-它们对应于信号的谐波。该功能将帮助我们做到这一点。 findpeaks 来自EngeeDSP函数集。让我们突出显示前六个谐波-这些峰值至少为-50dB，严重程度至少为30-并将其显示在图表上。:

spectrum, fvec = fftspectrum(sig, fs; nfft = nFFT, onesided = true);
pks, locs = EngeeDSP.Functions.findpeaks(spectrum, out=:data, 
                                                    MinPeakHeight = -50,
                                                    MinPeakProminence = 30);
plot(fvec, spectrum, linewidth=2, xlim=(0,1500),
        xguide = "Частота (Гц)", yguide = "Мощность (дБ)",
        title = "Гармоники на спектре аудиосигнала")
scatter!(fvec[locs], pks, legend = false)

第一个峰值为219Hz。但我们也可以通过编程方式确定前六个谐波的频率，并以Hz为单位找到它们之间的平均距离。我们将把这个估计值作为基频.:

harmonics = fvec[locs[1:6]]

6-element Vector{Float32}:
  219.0
  441.0
  661.0
  881.0
 1101.0
 1321.0

基频以统计方式计算为相邻峰值频率之差的算术平均值。:

base_freq = EngeeDSP.Functions.mean(diff(harmonics))

220.4f0

根据表,演奏的音符是一个小八度拉.

结论

我们熟悉了用于计算快速傅里叶变换的EngeeDSP功能，并实现了基于它估计信号频谱的算法，并解决了使用频谱分析确定音频信号特征的实际问题。