Основы цифровой обработки сигналов
07. Спектральный анализ
Спектральный анализ - одна из фундаментальных задач, лежащая в основе множества технологий: от сжатия аудио и изображений до радиолокации и медицинской диагностики.
Одна из основных целей спектрального анализа - оценка параметров сигнала, а именно:
- частоты
- амплитуды
- фазы
его гармонических составляющих.
В данном скрипте мы узнаем на примере функции fft из библиотеки EngeeDSP, как можно строить спектр и оценивать параметры тестового сигнала.
Простой сигнал и его спектр
Познакомимся с особенностями применения функции fft для анализа спектра на примере действительного косинусоидального сигнала с параметрами:
- длительность - 0.2 сек
- частота дискретизации - 2 кГц
- амплитуда - 0.5 В
- частота повторений -100 Гц
- начальная фаза - 60 градусов
fs = 2000;
dt = 1/fs;
t = 0:dt:0.2-dt;
pure_sine = 0.5 * cos.(2*pi*t*100 .+ pi/3);
plot(t,pure_sine,
leg=false,
title = "Чистая синусоида 100 Гц",
xguide = "Время (с)",
yguide = "Амплитуда (В)")
Стоит отметить длину нашего дискретного сигнала pure_sine- 400 отсчётов. Это важно для дальнейшего анализа.
nsamp = length(t)
Применим к нему функцию fft и визуализируем результат:
A = EngeeDSP.Functions.fft(pure_sine);
plot(A, leg=false)
Выход функции БПФ - это комплексный вектор A с таким же количеством отсчётов, как и у исходного сигнала во временной области.
Построение амплитудного спектра известного сигнала
Для оценки амплитудного спектра сигнала необходимо взять отсчёты выхода БПФ по модулю:
plot(abs.(A), leg=false)
Мы наблюдаем два пика, соответствующих положительной и отрицательной области частот. Но такое представление не даёт нам достоверной информации ни о частоте, ни об амплитуде этих пиков.
Для более явного отображения спектра воспользуемся функцией fftshift для сдвига выходного вектора БПФ, а также зададим вектор частот, который отложим по оси Х. Мы наблюдаем спектр сигнала в так называемой первой зоне Найквиста, которая простирается от -fs/2 до +fs/2. Разрешение по частоте (RBW, resolution bandwidth), или шаг частотной сетки, мы можем узнать, поделив весь рассматриваемый частотный диапазон (либо от -fs/2 до +fs/2, либо от 0 до fs) на количество точек БПФ:
B = EngeeDSP.Functions.fftshift(A);
println("Разрешение по частоте (RBW) = ", fs/nsamp, " Гц")
df = fs/nsamp;
freq_vec = -fs/2:df:fs/2-df;
plot(freq_vec, abs.(B),
xguide = "Частота (Гц)")
Теперь мы видим два пика на частотах -100 Гц и +100 Гц. Но их амплитуды не соответствуют амплитуде рассматриваемой синусоиды.
Для точной оценки энергетики сигнала на его спектре мы нормируем выходной вектор БПФ по количеству отсчётов, то есть делим на nsamp. Но так как половина энергии нашего действительного сигнала "ушла" в область отрицательных частот, нам надо умножить результирующий спектр на 2. Ну и для отображения воспользуемся линейчатым графиком (stem) с круглыми маркерами в области положительных частот от 0 до 400 Гц:
S = 2*B/nsamp;
println("Разрешение по частоте (RBW) = ", fs/nsamp, " Гц")
plot(freq_vec, abs.(S),
line=:stem,
marker=:circle,
xlim = (0, 400),
leg=false,
title = "Амплитудный спектр синусоиды",
xguide = "Частота (Гц)",
yguide = "Амплитуда (В)")
Мы наблюдаем пик на частоте 100 Гц амплитудой 0.5, что точно соответствует параметрам исходного синусоидального сигнала.
Фазовый спектр
Для оценки фазы колебания необходимо посчитать "угол" комплексного отсчёта функцией angle, и для удобства представления перевести его из радиан в градусы:
plot(freq_vec, angle.(S) * 180/pi,
line=:stem,
marker=:circle,
xlim = (0, 400),
leg=false,
title = "Фазовый спектр синусоиды",
xguide = "Частота (Гц)",
yguide = "Фаза (град)")
Изменение разрешения по частоте
Попробуем изменить количество отсчётов входного сигнала ("окно" БПФ). Подадим на вход функции fft первые 80 отсчётов синусоиды. В этом случае на выходе алгоритма БПФ так же получится 80 комплексных отсчётов, а значит изменится шаг частотной сетки. Сигнал стал короче в 5 раз, следовательно, разрешение по частоте стало 25 Гц:
winsize = 80;
rbw80 = fs/winsize;
freq80 = 0:rbw80:fs - rbw80;
S80 = (2*EngeeDSP.Functions.fft(pure_sine[1:winsize]))./winsize;
println("Разрешение по частоте (RBW) = ", rbw80, " Гц")
plot(freq80, abs.(S80),
line=:stem,
marker=:circle,
xlim = (0, 400),
title = "Амплитудный спектр синусоиды",
xguide = "Частота (Гц)",
yguide = "Амплитуда (В)")
Изменение частотной сетки на оценку амплитуды и частоты сигнала не повлияло. Дискрет в 25 Гц позволяет нам "точно попасть" в спектральный компонент на 100 Гц. Но что произойдёт, если шаг спектра не будет кратен частоте оцениваемой синусоиды?
winsize = 110;
rbw110 = fs/winsize;
freq110 = 0:rbw110:fs - rbw110;
S110 = (2*EngeeDSP.Functions.fft(pure_sine[1:winsize]))./winsize;
println("Разрешение по частоте (RBW) = ", rbw110, " Гц")
plot(freq110, abs.(S110),
line=:stem,
marker=:circle,
xlim = (0, 400),
title = "Амплитудный спектр синусоиды",
xguide = "Частота (Гц)",
yguide = "Амплитуда (В)")
Мы наблюдаем два пика на частотах 91 Гц и 109 Гц. Это явление утечки спектра (spectral leakage): энергетика сигнала на 100 Гц "расплылась" на соседние дискреты частотной сетки, между которыми она пришлась (более по-научному, была захвачена основным и боковыми лепестками оконной функции).
Дополнение нулями при БПФ
Zero padding (добавление нулей) при БПФ — это техника, которая заключается в дополнении исходного сигнала нулями до нужной длины перед выполнением БПФ. Добавление нулей увеличивает кол-во точек БПФ, уменьшая шаг частотной сетки и делая спектр визуально более гладким. Но это не добавляет новой информации о сигнале — просто интерполирует существующий спектр.
Дополним нашу синусоиду нулями до длины вектора в 2000 отсчётов:
winsize = 2000;
rbw2000 = fs/winsize;
freq2000 = 0:rbw2000:fs - rbw2000;
numofzeros = winsize - length(pure_sine);
println("Добавлено нулей = ", numofzeros)
S2000 = (2*EngeeDSP.Functions.fft([pure_sine; zeros(numofzeros)]))./length(pure_sine);
println("Разрешение по частоте (RBW) = ", rbw2000, " Гц")
plot(freq2000, abs.(S2000),
line=:stem,
marker=:circle,
xlim = (0, 200), leg=false,
title = "Амплитудный спектр синусоиды",
xguide = "Частота (Гц)",
yguide = "Амплитуда (В)")
Можно заметить, что при нормировке на размер выходного "окна" пик спектра находится на уровне 0.1, а не 0.5. Это связано с тем, что сигнал стал длиннее в пять раз за счёт нулей, и его энергетика, соответственно, также была ослаблена в пять раз. А значит надо нормировать на количество отсчётов реального ненулевого сигнала во времени.
Классический алгоритм работает эффективнее, когда число точек БПФ является степенью 2, и добавление нулей полезно для выравнивания длины и визуализации, но оно не увеличивает истинное разрешение по частоте. Для реального улучшения спектрального анализа нужно увеличивать длительность записи сигнала.
Интерактивный анализ спектра
В заключении, применим техники оценки спектра для сигнала, представляющего собой сумму трёх синусоид с частотами:
- 80 Гц
- 120 Гц
- 245 Гц
и амплитудами:
- 10 В
- 5 В
- 2.4 В
соответственно.
Отобразим форму результирующего вектора во временной области:
three_sines = [10 5 2.4] .* sin.(2*pi*t*[80 120 245]);
sig = sum(three_sines, dims = 2);
sig = vec(sig);
plot(t,sig,
leg=false,
title = "Сумма синусоид 80 Гц, 120 Гц и 245 Гц",
xguide = "Время (с)",
yguide = "Амплитуда (В)")
А теперь воспользуемся функционалом скриптов Engee для создания маски кодовой ячейки. Мы можем изменять значение разрешения по частоте (RBW) при помощи инструмента slider.
В коде ячейки происходит расчёт количества выходных точек БПФ, сравнение его с количеством точек дискретного сигнала во времени, и в случае превышения одного над другим выполняется или же дополнение нулями, или же усечение входного "окна" алгоритма. Затем вычисляется и отрисовывается амплитудный спектр.
RBW = 10 # @param {type:"slider",min:1,max:20,step:1}
Шум = 0 # @param {type:"slider",min:0,max:40,step:1}
nfft = Int64(round(fs/RBW));
println("Точек БПФ: ", nfft,
" Отсчётов сигнала: ", nsamp);
three_sines = [10 5 2.4] .* sin.(2*pi*t*[80 120 245]);
sumsines = vec(sum(three_sines, dims = 2));
sig = sumsines .+ Шум * randn(length(sig));
if nsamp >= nfft
win = sig[1:nfft];
Y = 2*EngeeDSP.Functions.fft(win)./nfft;
println("Усечение входного окна")
else
win = [sig; zeros(nfft - nsamp)];
Y = 2*EngeeDSP.Functions.fft(win)./nsamp;
println("Дополнение нулями: интерполяция спектра")
end
ff = (0:nfft-1)*fs/nfft;
plot(ff, abs.(Y),
line=:stem,
xlim = (0, 400),
marker =:circle,
lab = "Сумма синусоид",
title = "Амплитудный спектр сигнала",
xguide = "Частота (Гц)",
yguide = "Амплитуда (В)")
Попробуйте поэкспериментировать с различными значениями RBW, от 1 до 20 Гц, и посмотреть на точность оценки частоты и амплитуды отдельных спектральных компонент сигнала.
Вывод
Мы познакомились с особенностями применения функции fft библиотеки EngeeDSP для оценки спектра дискретных сигналов, а также узнали, как применять макси кодовых ячеек и выполнять задачи цифровой обработки сигналов в интерактивном режиме.