Основы цифровой обработки сигналов
Спектральное представление цифровых сигналов
Частотное представление сигнала
Частотное представление сигнала – это зависимость измеряемых параметров сигнала от частоты.
Зависимость энергии сигнала от частоты называется спектром сигнала.
Спектр сигнала может быть непрерывным (как показано на верхнем рисунке) или дискретным (как показано на нижнем рисунке).
В основе частотного представления сигнала лежит разложение сложного сигнала на сумму простых сигналов. Сложный сигнал представляется в виде суммы синусоид с разными частотами, умноженных на весовые коэффициенты. Весовые коэффициенты удобно представлять комплексными числами: модуль числа соответствует амплитуде, а аргумент – начальной фазе.
Действительная синусоида в частотной области
Рассмотрим математическое описание дейтвительной синусоиды
Разложим ее по формуле Эйлера:
или
Мы видим два слагаемых с зеркальными начальными фазами и частотами. Следовательно, на частотной оси мы получим две гармоники с одинаковыми амплитудами, одинаковыми по значению, но разными по знаку частотами и начальными фазами:
Спектр действительной синусоиды выглядит так:
Спектр любого действительного сигнала всегда симметричен относительно нулевой частоты.
Преобразование Фурье
Преобразование Фурье используется как способ разложения сигнала на частоты и амплитуды, т.е. для перевода сигнала из временного представления в частотное.
Для непрерывных сигналов применяется интегральная форма преобразования Фурье (непрерывное преобразование Фурье):
а для дискретных сигналов – дискретное преобразование Фурье (ДПФ):
Результатом ДПФ будет набор комплексных отсчетов спектра.
Существуют обратные преобразования Фурье (непрерывное и дискретное), переводящие сигнал из частотной области во временную.
Быстрое преобразование Фурье (БПФ) – это эффективный алгоритм вычисления ДПФ. Оптимальная длина временной последовательности для БПФ равна , где – натуральное число.
Наиболее распространенная версия БПФ – алгоритм Кули–Тьюки. Он основывается на том факте, что для вычисления ДПФ из, например, 8 точек, оно может быть выполнено отбором четных и нечетных значений входного сигнала и выполнением двух четырехточечных ДПФ с дальнейшей операцией прореживания во времени.
Пример. Создадим в Engee аккорд ля мажор первой октавы, состоящий из гармонических колебаний с частотами 440, 550 и 660 Гц.
using Plots, FFTW, DSP;
fc = [440 550 660]';
fs = 8000;
dt = 1/fs;
t = [0:dt:0.1;];
x = cos.(2*pi*fc[1]*t) + cos.(2*pi*fc[2]*t) + cos.(2*pi*fc[3]*t);
x = x/3;
Выполним преобразование Фурье функцией fft, которая содержится в библиотеке FFTW.
X = fft(x);
Построим амплитудный спектр сигнала. Воспользуемся функцией abs для вычисления модуля комплексного вектора X. Для вычисления вектора частот воспользуемся функцией fftfreq, аргументами которой являются число отсчетов сигнала и частота дискретизации.
X = abs.(X);
fr = fftfreq(length(t), fs);
plot(fr, X, xlabel="частота, Гц", ylabel="амплитудный спектр", legend=false)
Применение БПФ
БПФ и обратное БПФ применяются во многих задачах ЦОС:
- Фильтрация в частотной области;
- Сжатие;
- "Быстрая" дискретная свертка;
- Спектральный анализ.
Спектральный анализ сигнала
При выполнении БПФ выходной вектор называется точками, или отсчетами БПФ. Отсчеты сигнала обычно поступают скалярно (по одному), поэтому перед БПФ осуществляется буферизация – накопление отсчетов вектора. Входной вектор временных отсчетов называется окном. Размер входного окна тем или иным способом приравнивается к длине БПФ (). Если во входном окне отсчетов больше, чем нужно, то лишние отбрасываются. Если меньше, то входной вектор дополняется нулями.
Разрешение по частоте зависит от длины БПФ и от частоты дискретизации сигнала:
Берется весь частотный диапазон от до и заполняется точками, в которых нужно оценить спектр. Расстояние между соседними точками – это и есть разрешение по частоте.
Оконные функции при спектральном анализе
Спектральная плотность мощности – мощность сигнала, приходящаяся на единичный интервал частоты.
Часто мы не можем определить точное положение той или иной гармоники в спектре, т.к. она попадает между точками БПФ. Энергия этой гармоники может проявиться в соседних точках БПФ. Этот эффект называется утечкой спектра. Избежать его полностью невозможно, но сгладить его влияние можно с помощью оконных функций.
Оконная функция – это набор весовых коэффициентов, изменяющих значения амплитуды отсчетов сигнала во временной области.
При анализе большого сигнала мы разбиваем его на небольшие отрезки, которые мы взвешиваем оконной функцией и отправляем на вход алгоритма БПФ.
Выполним в Engee спектральный анализ при помощи встроенных функций. Входные данные остались теми же (аккорд ля мажор первой октавы). Воспользуемся встроенной функцией welch_pgram для оценки спектральной плотности мощности методом Уэлча. Входными аргументами этой функции являются размер входного окна, длина БПФ, количество перекрывающихся отсчетов и тип оконной функции. В нашем случае выбрана функция Хэмминга.
n=div(length(x),8);
y = welch_pgram(x, n, div(n,2), onesided=true, nfft=length(x), fs=fs, window=DSP.hamming);
plot(freq(y), power(y), xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
✏️Задание 1
Сформируйте в Engee сигнал, представляющий собой синусоиду на частоте 1000 Гц с добавлением аддитивного белого гауссовского шума ( дБ). Постройте спектр сигнала двумя способами: с помощью функции fft и welch_pgram.
Решение
using Plots, DigitalComm, FFTW;
fs = 4000;
t = [0:1/fs:0.2;];
sin_wave = cos.(2*pi*1000*t);
x = addNoise(sin_wave, 5)[1];
X = abs.(fft(x));
fr = fftfreq(length(t), fs);
plot(fr, X, xlabel="частота, Гц", ylabel="амплитудный спектр", legend=false)
n=div(length(x),8);
y = welch_pgram(x, n, div(n,2), onesided=true, nfft=length(x), fs=fs, window=DSP.hamming);
plot(freq(y), power(y), xlabel="частота, Гц", ylabel="спектральная плотность мощности", legend=false)
Нелинейные искажения
Линейными называются изменения сигнала, приводящие к изменению амплитуд и/или фаз его спектральных компонент.
Нелинейными называются изменения сигнала, приводящие к появлению новых спектральных компонент.
Оконное преобразование Фурье
Оконное преобразование Фурье (short-time Fourier transform, STFT) – это метод оценки сигнала, спектр которого изменяется во времени.
Оконное БПФ позволяет получать кусочные спектры, которые затем накладываются на временную ось.
Спектрограммой называется трехмерная картина изменения спектра во времени.
Пример. Построим в Engee спектрограмму сигнала, представляющего собой синусоиду на частоте 1000 Гц c амплитудой, возрастающей с течением времени по закону . Для этого воспользуемся встроенной функцией spectrogram, содержащейся в библиотеке DSP. Входные аргументы этой функции полностью аналогичны входным аргументам функции welch_pgram. Выход функции spectrogram представляет собой объект с тремя полями:
time– вектор времени,freq– вектор частот,power– матрица спектральной плотности мощности.
С помощью функции plot отобразим спектрограмму на графике. Величина спектральной плотности мощности переведена в децибелы функцией pow2db. На графике она показана различной интенсивностью цвета.
using Plots, DSP;
fs = 4000;
t = [0:1/fs:0.2;];
x = (t.^3).*cos.(2*pi*1000*t);
n=div(length(x),8);
y = DSP.spectrogram(x, n, div(n,2), onesided=true, nfft=length(x), fs=fs, window=DSP.hamming);
plot(y.time, y.freq, pow2db.(y.power), xlabel="частота, Гц", ylabel="t, с", legend=false)
✏️Задание 2
Сформируйте в Engee сигнал из Задания 1 – синусоиду на частоте 1000 Гц с добавлением аддитивного белого гауссовского шума ( дБ). Постройте и выведите на экран спектрограмму этого сигнала.
Решение
using Plots, DigitalComm, DSP;
fs = 4000;
t = [0:1/fs:0.2;];
sin_wave = cos.(2*pi*1000*t);
x = addNoise(sin_wave, 5)[1];
n=div(length(x),8);
y = DSP.spectrogram(x, n, div(n,2), onesided=true, nfft=length(x), fs=fs, window=DSP.hamming);
plot(y.time, y.freq, pow2db.(y.power), xlabel="частота, Гц", ylabel="t", legend=false)