День 3 Летней школы Julia
Часть 3: спектральный анализ
В цифровой обработке сигналов спектральный анализ — это набор методов, позволяющих разложить сигнал на частотные компоненты и исследовать его свойства в частотной области.
Мы рассмотрим методы оценки спектра на основе БПФ, периодограммы и Уэлча для задачи анализа аудиосигнала.
# Это нужно выполнить, если не установлены необходимые библиотеки
#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
Знакомство с БПФ
В качестве первого тестового сигнала возьмём сумму трёх синусоид из первой части:
fs = 1000; # частота дискретизации сигнала
dt = 1/fs; # период дискретизации
stoptime = 1; # время окончания сигнала
t = 0:dt:stoptime-dt; # вектор отсчётов времени
three_sines = [1.4 0.6 1] .* sin.(2*pi*t.*[10 15 30] .+ [pi/4 0 pi/6]);
sum_sines = vec(sum(three_sines, dims = 2));
plot(t, sum_sines, title = "Сумма трёх синусоид")
Попробуем применить к сигналу функцию fft и визуализировать результат:
A = fft(sum_sines);
plot(A)
Комплексный вектор на выходе БПФ содержит информацию как об амплитуде составляющих, так и о фазе. Отобразим выход БПФ по модулю комплексного числа:
plot(abs.(A))
Функция fftshift позволяет представить выход БПФ в привычном виде спектра, в нашем случае симметричного относительно центра - нулевой частоты. Но нам также придётся отложить по оси абсцисс вектор дискретных частот в пределах от -fs/2 до fs/2:
B = fftshift(A);
nsamp = length(t);
df = fs/nsamp;
println("Разрешение по частоте = ", df, " Гц")
freq_vec = -fs/2:df:fs/2 - df;
plot(freq_vec, abs.(B))
Наконец, чтобы получить амплитудный спектр сигнала, нам необходимо учитывать количество отсчётов на выходе функции БПФ, а также принять во внимание, что энергетика вещественного сигнала также распространяется и на область отрицательных частот, которая, в нашем случае, неинформативна:
S = 2*B/nsamp;
plot(freq_vec, abs.(S), line=:stem, marker=:circle, xlim=(0, 100))
Исходный сигнал для анализа
В качестве тестового сигнала мы берём аудиозапись гитарной струны. Наша задача - определить, что за нота проигрывается, анализируя спектр сигнала.
Загрузим численные отсчёты сигнала функцией wavread, которая также позволяет выделить частоту дискретизации сигнала fs.
audio_data, fs = wavread("guitar.wav");
Подсчитаем количество отсчётов сигнала во временной области, вычислим временной шаг и сформируем вектор времени:
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)
Узнать основную частоту колебания по форме периодического сигнала во временной области можно, но всё же весьма затруднительно.
Оценка спектра сигнала
Рассмотрим три метода для построения спектра мощности или графика спектральной плотности мощности (СПМ) тестового сигнала. Для качественного анализа частотного состава подходят обе метрики.
Начнём с построения спектра мощности на основе БПФ:
A = fft(sig);
B = fftshift(A);
S = (2*B)/nsamples;
df = fs/nsamples;
freq_vec = -fs/2:df:fs/2-df;
plot(freq_vec, pow2db.(abs.(S.^2)))
На спектре явно прослеживается гармоническая структура - равнорасположенные пики. Гармоники - это синусоидальные колебания, частоты которых кратны основной частоте .
- самая низкая частота в сигнале (первая гармоника), а высшие гармоники - частоты вида , где (вторая, третья и т. д.).
Периодограмма
Периодограмма — это метод оценки спектральной плотности мощности (СПМ) сигнала на основе его дискретного преобразования Фурье (ДПФ). На практике вычисляется через БПФ.
Алгоритм:
-
Берётся отрезок сигнала (например, 1024 отсчёта).
-
Вычисляется БПФ этого отрезка.
-
Амплитудный спектр возводится в квадрат и нормируется на длину выборки.
-
Получается график мощности по частотам.
p1 = DSP.periodogram(sig, onesided = true, fs = fs);
plot(freq(p1), pow2db.(power(p1)))
Метод Уэлча
Метод Уэлча - это модификация периодограммы, в которой:
-
Сигнал разбивается на перекрывающиеся отрезки.
-
К каждому отрезку применяется оконная функция (например, Ханна, Хэмминга) - это уменьшает эффект утечки спектра.
-
Для каждого отрезка считается периодограмма.
-
Результаты усредняются → снижается дисперсия.
p2 = DSP.welch_pgram(sig, onesided = true, fs = fs);
plot(freq(p2), pow2db.(power(p2)))
Определение основной частоты сигнала
Остановимся на оценке спектра методом Уэлча и проанализируем область низких частот (от 0 до 1500 Гц). Найдём пики в спектре - они соответствуют гармоникам сигнала. В этом нам поможет функция findpeaks1d из соответствующего пакета. Выделим первые шесть гармоник, и отобразим их на графике:
spectrum = pow2db.(power(p2));
fvec = freq(p2);
idx, properties = findpeaks1d(spectrum; height = -65, prominence = 20);
plot(fvec, spectrum, linewidth=3, xlim=(0,1500))
scatter!(fvec[idx], spectrum[idx])
Самый первый пик находится на частоте 221 Гц. Но мы также можем программно определить частоты первых шести гармоник и найти среднее расстояние между ними в Гц. Именно эту оценку мы и примем за основную частоту:
harmonics = fvec[idx[1:6]]
base_freq = mean(diff(harmonics))
Судя по таблице, сыгранная нота - это Ля малой октавы.