DSP基础：频谱分析
在**数字信号处理(DSP)*频谱分析是一组方法，使得能够将信号分解成频率分量并研究其在频域中的属性。
DSP中频谱分析的主要任务:
-
信号的频率组成的确定（例如，音频信号中音调分量的分配）。
-
隐藏周期的检测（例如，机械系统中的振动分析）。
-
滤波器规格（用于随后的噪声去除和有用的频率提取）。
-
数据压缩（使用JPEG，MP3等频率表示）。
-
模式识别（例如，语音信号分析）。
分析的主要方法有:
*离散傅立叶变换（dft）和快速傅立叶变换（FFT）。 了解更多关于FFT-在这个примере.
*周期图(功率谱估计)
*韦尔奇法（平均周期图）
*参数化方法(autoregression,AR)
*过滤器库
在本例中，我们将考虑用于音频信号分析的FFT、周期图和Welch方法。
# Это нужно выполнить, если не установлены необходимые библиотеки
#Pkg.add("DSP")
#Pkg.add("WAV")
#Pkg.add("Base64")
#Pkg.add("FFTW")
#Pkg.add("FindPeaks1D")
#Pkg.add("Statistics")
using DSP, WAV, Base64, FFTW, FindPeaks1D, Statistics
分析的初始信号
作为测试信号，我们使用吉他弦的录音。 我们的任务是通过分析信号的频谱来确定正在播放哪个音符。
让我们用函数加载信号的数值样本 wavread，这也允许您选择信号的采样频率 fs.
filename = "$(@__DIR__)/guitar.wav"
audio_data, fs = wavread(filename);
让我们统计信号在时域上的样本数，计算时间步长并形成时间向量。:
nsamples = length(audio_data);
dt = 1/fs;
t = 0:dt:((nsamples - 1) .* dt);
sig = vec(audio_data);
我们将使用一个附加功能来收听音频。 audioplayer:
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(sig,fs)
最后，让我们在时域中显示信号。:
plot(t, sig,
xguide = "Время, с",
yguide = "Амлпитуда, В",
title = "Аудиосигнал во времени",
legend = false)
在时域上通过周期信号的形状找出基波振荡频率是可能的，但仍然非常困难。
信号频谱的估计
让我们考虑三种方法来构建功率谱或测试信号的频谱功率密度（SPM）图。 这两种度量都适用于频率组成的定性分析。
让我们从基于FFT构建功率谱开始。 了解有关函数和数学的更多信息-在前面的примере.
A = fft(sig);
B = fftshift(A);
df = fs/nsamples;
freq_vec = -fs/2:df:fs/2-df;
plot(freq_vec, pow2db.(abs.(B.^2)),
xguide = "Частота, Гц",
yguide = "Мощность, дБ",
title = "Спектр мощности сигнала (на основе БПФ)",
xlim = (0,10000),
legend = false)
在光谱上清晰可见谐波结构-等间距的峰。 谐波是正弦振荡，其频率是基频的倍数。 .
-在信号的最低频率（一次谐波），和高次谐波是频率的形式 ，在哪里 （第二，第三等。).
周期图
Periodogram是一种用于估计频谱功率密度(SPM)基于其离散傅里叶变换(dft)的信号的方法。 在实践中，它是使用FFT计算的。
算法:
-
取信号的一段（例如，1024个样本）。
-
计算该段的FFT。
-
幅度谱被平方并由样本长度归一化。
-
得到频率功率图。
p1 = DSP.periodogram(sig, onesided=true, fs=fs);
plot(freq(p1), pow2db.(power(p1)),
xguide = "Частота, Гц",
yguide = "Мощность/частота, дБ/Гц",
title = "Спектральная плотность мощности сигнала (периодограмма)",
xlim = (0,10000),
legend = false)
然而，这种方法有许多缺点。:
*高方差-估计值随着小数据变化而"跳跃"。
*频谱的泄漏（泄漏）–由于信号的有限长度，会发生虚假频率（通过窗口函数解决，例如，Hanna或Hamming）。
*低频率分辨率-如果信号短，附近的频率合并。
韦尔奇方法
韦尔奇方法是对周期图的修改，其中:
*信号被分割成重叠的段。
*窗口函数应用于每个段（例如，Hanna，Hamming）-这减少了频谱泄漏的影响。
*为每个段计算周期图。
*结果平均→方差减少。
p2 = DSP.welch_pgram(sig, onesided=true, fs=fs);
plot(freq(p2), pow2db.(power(p2)),
xguide = "Частота, Гц",
yguide = "Мощность/частота, дБ/Гц",
title = "Спектральная плотнотсть мощности сигнала (метод Уэлча)",
xlim = (0,10000),
legend = false)
确定信号的基频
让我们专注于使用Welch方法评估SPM，并分析低频范围（从0到1500Hz）。 我们将在频谱中找到峰值-它们对应于信号的谐波。 该功能将帮助我们做到这一点。 findpeaks1d 从相应的包。 让我们选择前六个谐波并将它们显示在图表上。:
spectrum = pow2db.(power(p2));
fvec = freq(p2);
pkindices, properties = findpeaks1d(spectrum;
height = -65,
prominence = 20);
plot(fvec, spectrum,
linewidth = 3,
xguide = "Частота, Гц",
yguide = "Мощность/частота, дБ/Гц",
title = "Пики спектра сигнала (метод Уэлча)",
xlim = (0,1500),
label = "спектр")
scatter!(fvec[pkindices], spectrum[pkindices]; color="red", markersize=5, label="пики")
第一个峰值为219Hz。 但我们也可以通过编程方式确定前六个谐波的频率，并以Hz为单位找到它们之间的平均距离。 我们将把这个估计值作为基频.:
harmonics = fvec[pkindices[1:6]]
base_freq = mean(diff(harmonics))
根据表,演奏的音符是一个小八度拉.
选择第一谐波
我们合成了一个简单的椭圆低通滤波器（低通滤波器），可以抑制300Hz以上的频率。 这将使我们能够切断高次谐波。 让我们过滤掉原始信号并听取结果。:
myfilt = digitalfilter(Lowpass(2*300/fs), Elliptic(10, 0.1, 80));
clean220 = filt(myfilt, sig);
audioplayer(clean220 * 2, fs)
我们听到更"纯粹"的音调，没有色彩。 不能立即识别其中的吉他，因为它是高次谐波给乐器音色。
结论
我们研究了三种估计信号频谱的标准方法-FFT，周期图和Welch方法，使用分析来自振荡吉他弦的音频信号的频率组成的例子-并且还确定了振荡的基频。