#include <stdio.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
#define N 100
#define P 10000
#define F0 100.0
#define F1 385.0
#define FS 1.0E3
#define STR_LEN 1024
/*******************************************************************************
Расчет спектральной плотности мощности сигнала, содрежащего 2 комплексные
гармоники на частотах F0 и F1
s(n) = exp(j * 2 * pi * n * F0 / Fs) + exp(j * 2 * pi * n * F1 / Fs) + r(n),
n = 0 ... N-1; r(n) - белый гауссов шум
*******************************************************************************/
void scalloping(int argc, char* argv[], int win_type, char* fn)
{
complex_t s[N], y[N];
double f[N], S[N], A[N], t[N], w[N], wn;
int n, k;
char fname[STR_LEN] = {0};
char pngname[STR_LEN] = {0};
char cmd[STR_LEN] = {0};
fft_t pfft = {0}; /* FFT structre */
random_t rnd = {0}; /* Random structur */
void* hplot; /* GNUPLOT handle */
/* Инициализация генератора белого шума */
random_init(&rnd, RAND_TYPE_MT19937, NULL);
/* время t = n = 0 ... N-1 */
linspace(0, N, N, DSPL_PERIODIC, t);
/* генерируем окно */
window(w, N, win_type, 0);
/* расчет wn = W(0) = sum[w(n)] */
sum(w, N, &wn);
/* формируем сигнал содержащий 2 комплексные
экспоненты на частотах F0 и F1 */
for(k = 0; k < N; k++)
{
RE(s[k]) = cos(M_2PI * t[k] * F0 / FS);
IM(s[k]) = sin(M_2PI * t[k] * F0 / FS);
RE(s[k]) += cos(M_2PI * t[k] * F1 / FS);
IM(s[k]) += sin(M_2PI * t[k] * F1 / FS);
}
/* Массив А - оценка спектральной плотности мощности рассчитывается как
усредненная по P реализациям оценка спектральной плотности энергии */
memset(A, 0, N*sizeof(double));
/* Усреднение спектральных плотностей энергии */
for(n = 0; n < P; n++)
{
/* формируем реализацию комплексного белого шума на уровне 1 дБ / Гц */
randn((double*)y, 2*N, 0.0, 1.122/sqrt(2), &rnd);
/* Добавляем к сигнлу реализацию шума и умножаем на оконную функцию:
y(n) = w(n) * [s(n) + r(n)] */
for(k = 0; k < N; k++)
{
RE(y[k]) = (RE(y[k]) + RE(s[k])) * w[k] * N / wn;
IM(y[k]) = (IM(y[k]) + IM(s[k])) * w[k] * N / wn;
}
/* Рассчитываем спектральную плотность энергии с перестановкой
результата ДПФ симметрично относительно 0 частоты */
fft_mag_cmplx(y, N, &pfft, FS, DSPL_FLAG_FFT_SHIFT, S, f);
/* Добавляем текущую оценку к массиву A */
for(k = 0; k < N; k++)
A[k] += S[k];
}
/* Усредняем накопленный массив A, нормируем к размеру N
и переводим в логарифмическую шкалу */
for(k = 0; k < N; k++)
A[k] = 10.0*log10(A[k] /(P*N));
/* сохраняем в файл dat/sclloping_[окно].txt */
memset(fname, 0, STR_LEN);
sprintf(fname, "dat/sclloping_%s.txt", fn);
writetxt(f, A, N, fname);
/* Очищаем структуру FFT */
fft_free(&pfft);
/* plotting by GNUPLOT */
/* Create window plot */
memset(pngname, 0, STR_LEN);
sprintf(pngname, "dat/sclloping_%s.png\n", fn);
gnuplot_create(argc, argv, 420, 420, pngname, &hplot);
gnuplot_cmd(hplot, "set xlabel 'f'");
gnuplot_cmd(hplot, "set ylabel '|S(f)|^2, dB'");
gnuplot_cmd(hplot, "unset key");
gnuplot_cmd(hplot, "set xrange[50:450]");
gnuplot_cmd(hplot, "set yrange[-2:22]");
memset(cmd, 0, STR_LEN);
sprintf(cmd, "set title '%s'\n", fn);
gnuplot_cmd(hplot, cmd);
memset(cmd, 0, STR_LEN);
sprintf(cmd, "plot \'%s\' with impulses\n", fname);
gnuplot_cmd(hplot, cmd);
gnuplot_close(hplot);
}
/*******************************************************************************
Основная программа
*******************************************************************************/
int main(int argc, char* argv[])
{
void* hdspl;
/* Загрузка DSPL-2.0 */
hdspl = dspl_load();
if(!hdspl)
{
printf("cannot to load libdspl!\n");
return 0;
}
/* Спектральная плотность мощности при использовании прямоугольного окна */
scalloping(argc, argv, DSPL_WIN_RECT | DSPL_WIN_PERIODIC, "rect");
/* Спектральная плотность мощности при использовании окна Хемминга */
scalloping(argc, argv, DSPL_WIN_HAMMING | DSPL_WIN_PERIODIC, "hamming");
/* Спектральная плотность мощности при исп. окна Блэкмана - Харриса */
scalloping(argc, argv, DSPL_WIN_BLACKMAN_HARRIS | DSPL_WIN_PERIODIC,
"blackmanharris");
/* free dspl handle */
dspl_free(hdspl);
return 0;
}
Функции оконного сглаживания
![]() DSPL-2.0 — свободная библиотека алгоритмов цифровой обработки сигналов Распространяется под лицензией LGPL v3
Страница проекта на SourceForge
|
В предыдущем разделе мы анализировали эффект растекания спектра дискретного преобразования Фурье (ДПФ) и говорили, что он возникает в результате дискретизации дискретно-временного преобразования Фурье (ДВПФ) ограниченного во времени сигнала.
Пара ДВПФ сигнала ,
имеет вид:
![equation 1](img/eq-01.png)
ДВПФ спектр дискретного сигнала представляет собой -периодическую функцию нормированной частоты
.
Взяв дискретных значение
на сетке частот
,
получим пару дискретного преобразования Фурье:
![equation 2](img/eq-02.png)
Рассмотрим
дискретно-временное преобразование Фурье
комплексного гармонического сигнала , ограниченной
длительности
отсчетов, частота
изменяется в пределах от
до
. Ограничение по длительности сигнала
осуществляется умножением на оконную функцию
. Тогда ДВПФ
равно:
![equation 3](img/eq-03.png)
Таким образом, ДВПФ комплексного гармонического сигнала полностью описывается ДВПФ оконной функции, которая используется для ограничения длительности сигнала.
В
предыдущем разделе
мы рассмотрели эффект растекания ДПФ, который возникает, когда частота гармонического сигнала не совпадает с частотами дискретизации ДВПФ. Мы отмечали, что эффект растекания спектра приводит к гребешковым искажениям амплитуды гармонического сигнала, а также к уменьшению динамического диапазона спектрального анализа ввиду маскирования слабых сигналов под боковыми лепестками мощных гармоник при растекании спектра последних.
Эффект растекания ДПФ возникает в результате дискретизации ДВПФ ограниченного во времени сигнала на боковых лепестках спектра оконной функции . Для того чтобы уменьшить растекание, надо уменьшить уровень боковых лепестков спектра
.
Как мы знаем, уровень боковых лепестков спектра и скорость их убывания по частоте зависит от степени гладкости исходной функции. Поэтому для уменьшения эффекта растекания необходимо изменить оконную функцию , исключив скачки в начале и конце. На рисунке показаны
в виде прямоугольного окна (несглаженная) и сглаженное окно Хэмминга, а также их спектральные плотности энергии
.
![ДВПФ прямоугольного окна и окна Хемминга,](img/win_ex0.png)
![N=50](img/eqlin-017.png)
Можно видеть, что уровень боковых лепестков спектра окна Хемминга значительно ниже чем у прямоугольного окна. Таким образом, использовав сглаженные окна, можно получить увеличение динамического диапазона спектрального анализа, и также уменьшить уровень гребешковых искажений.
В спектральном анализе спектры ДВПФ и ДПФ могут выводится на график и представляется в различных единицах частоты. Например ДВПФ определяется как
периодическая функция нормированной частоты
. Тогда один период ДВПФ может быть переставлен в диапазоне от 0 дo
. В этом случае, ввиду периодичности
, частотный диапазон от
до
соответствует отрицательной области от
до нуля, как это показано на рисунке 2а.
Если для вывода спектра на график необходимо разместить в центре, то можно представить частоту
от
до
и произвести перестановку спектральных компонент, как это показано на рисунке 2б.
![Пересчет единиц частоты](img/win_freq_recalc.png)
От нормированной частоты ДВПФ можно перейти к спектральному представлению
, где
— частота в единицах герц. Для этого необходимо знать частоту дискретизации исходного сигнала
, которая использовалась для оцифровки данных. Тогда
Гц, как это показано на рисунках 2б и г.
Ну и наконец, для анализа характеристик оконной функции мы будем нормировать частоту к частоте
в единицах спектральных отсчетов ДПФ (ДПФ бин). В этом случае нормировка задается выражением
ДПФ бин, где
— размер выборки данных (количество точек ДПФ). Тогда
ДПФ бин соответсвует первому отсчету ДПФ
. Графически шкала частот в единицах ДПФ бин показана на рисунках 2д и е.
Можно многими способами сгладить прямоугольное окно, и таким образом, получить множество семейств различных оконных функций. Поэтому нам надо рассмотреть параметры оконных функций, которые так или иначе влияют на качество спектрального анализа. Имея набор параметров мы сможем сравнить окна между собой и получить объективный подход для выбора того или иного окна для практического применения.
При совпадении частоты с частотой некоторого бина ДПФ имеет место соотношение
, где
— цело число от
до
. Тогда ДПФ
и имеет максимально-возможное значение.
В случае прямоугольного окна
при
, тогда
и имеет место максимальное когерентное усиление равное
среди всех возможных окон. В этом случае комплексная экспонента соответствующая
-му бину ДПФ представляет собой согласованный фильтр, для гармонического сигнала на частоте
и мы имеем максимальные отклик на выходе согласованного фильтра. Таким образом, величина когерентного усиления на выходе ДПФ для комплексной экспоненты соответствующей
-му бину ДПФ равно
раз по амплитуде или
![equation 4](img/eq-04.png)
![G_C = N](img/eqlin-035.png)
![|W(0)|/N](img/eqlin-036.png)
![equation 5](img/eq-05.png)
Еще раз заметим, что когерентное усиление относится только к случаю, когда частота гармонического сигнала попадает на бин ДПФ и эффекта растекания не наблюдается.
Если частота гармонического сигнала не совпадает с частотой бина ДПФ, то наблюдаются гребешковые искажения амплитуды гармонического сигнала на выходе ДПФ, ввиду того, что дискретизация ДВПФ производится на боковых лепестках спектра оконной функции
.
Поскольку ДВПФ оконной функции является симметричной функцией относительно нулевой частоты, то смещение частоты
на половину бина ДПФ приведет к максимальному уровню гребешковых искажений
, соответствующему частоте
.
Для исключения зависимости от длины выборки можно нормировать уровень гребешковых искажений к уровню когерентного усиления. Тогда
![equation 6](img/eq-06.png)
![equation 7](img/eq-07.png)
На рисунке 3 показаны частотные параметры спектральной плотности энергии оконных функций.
![Спектральные характеристики оконной функции](img/win_freq.png)
Значение частоты на рисунке 3 приведено в единицах спектральных отсчетов ДПФ (бин ДПФ).
На графиках можно видеть значение когерентного усиления (в данном случае усиление не нормировано к значению усиления прямоугольного окна). Относительно
можно выделить уровень максимального бокового лепестка
оконной функции.
Главный лепесток оконной функции характеризуется шириной по уровню дБ и
дБ (полосы
и
соответственно), а также шириной главного лепестка по нулевому уровню
.
Значение на частоте
бина ДПФ показывает максимальный уровень гребешковых искажений
относительно уровня
.
![W(\omega)](img/eqlin-014.png)
![\mathrm{w}(n)](img/eqlin-013.png)
![N](img/eqlin-005.png)
![equation 8](img/eq-08.png)
![\omega](img/eqlin-004.png)
![W(\omega)](img/eqlin-014.png)
![2\pi](img/eqlin-003.png)
![equation 9](img/eq-09.png)
![\mathrm{ENBW}](img/eqlin-049.png)
![\mathrm{ENBW}](img/eqlin-049.png)
![|W(\omega)|^2](img/eqlin-016.png)
![|W(0)|^2](img/eqlin-050.png)
![Эквивалентная шумовая полоса оконной функции](img/win_enbw.png)
![equation 10](img/eq-10.png)
![equation 11](img/eq-11.png)
![equation 12](img/eq-12.png)
![equation 13](img/eq-13.png)
Тогда можно заключить, что
![equation 14](img/eq-14.png)
![equation 15](img/eq-15.png)
![\mathrm{ENBW}](img/eqlin-049.png)
![\omega = 0](img/eqlin-019.png)
![\omega = 2\pi](img/eqlin-051.png)
![N](img/eqlin-005.png)
![equation 16](img/eq-16.png)
![\mathrm{w}(n) = 1](img/eqlin-032.png)
![n = 0\ldots N-1](img/eqlin-002.png)
![\mathrm{ENBW}= 1](img/eqlin-052.png)
![|W(\omega)|^2](img/eqlin-016.png)
![L_P](img/eqlin-053.png)
Поскольку эквивалентная шумовая полоса сглаженных окон выше чем прямоугольного окна, то это означает более высокий уровень шума на выходе.
Это можно пояснить так. Поскольку прямоугольное окно обладает максимальным когерентным усилением и является согласованным фильтром для гармонического сигнала на частоте бина ДПФ, то оно обеспечивает максимальное отношение сигнал–шум на выходе. Любое другое сглаженное окно будет обеспечивать худшее отношение сигнал–шум для гармонического сигнала на частоте бина ДПФ.
На рисунке 5 показан результат оценки спектральной плотности мощности сигнала , содержащего две гармоники на частотах
Гц и
Гц на фоне белого гауссова шума при частоте дискретизации 1 кГц и использовании ДПФ размера
точек и различных оконных функциях.
Оценка велась при усреднении по ансамблю из 10000 реализаций.
![Потери при обработке гармонических сигналов](img/scalloping_loss.png)
Заметим, что при построении графиков, разница в когерентном усилении оконных функций была скомпенсирована и выровнена с усилением прямоугольного окна. В результате значение гармоник на частоте 100 Гц для всех окон равно 20 дБ.
Можно видеть, что прямоугольное окно обеспечивает максимальное отношение сигнал–шум дБ для сигнала на частоте 100 Гц, но при этом наблюдаются максимальные гребешковые искажения
дБ.
При использовании окна Хемминга величина гребешковых искажений уменьшилась до значения дБ, однако отношение сигнал–шум также ухудшилось на величину
дБ.
Для окна Блэкмана-Харриса гребешковые искажения составляют дБ, в то время как отношение сигнал–шум ухудшается на величину
дБ.
Можно обратить внимание, что суммарные потери на обработку, которые включают гребешковые искажения и ухудшения отношения сигнал–шум для наихудшего случая, т.е. сигнала на частоте 385 Гц составляют дБ для прямоугольного окна,
дБ для окна Хемминга и
дБ для окна Блэкмана-Харриса. Таким образом, в наихудшем случае, суммарные потери на обработку при использовании прямоугольного окна, как правило, остаются выше, чем при использование сглаженных окон[1].
Данный пример был рассчитан с использованием библиотеки DSPL–2.0. Исходный код программы scalloping_loss.c расчета данных для построения графиков рисунка 5 приведён в следующем листинге.
Оконные функции используются в различных задачах спектрального анализа, а также при синтезе цифровых КИХ-фильтров. В задачах спектрального анализа мы стремимся получить более узкий основной лепесток спектра оконной функции, в то время как в задачах синтеза КИХ-фильтров часто требуется обеспечить линейность ФЧХ, что требует симметрии оконной функции во временно́й области.
Для удовлетворения требований спектрального анализа и задач синтеза линейнофазовых КИХ-фильтров различают периодические и симметричные оконные функции. Периодические оконные функции обладают лучшими частотными свойствами, но не имеют свойства временно́й симметрии. Пример периодических и симметричных треугольных окон Бартлетта показан на рисунке 6 для четной и нечетной
длин окна.
![Периодические и симметричные окна Бартлетта для и](img/win_periodic.png)
![N = 8](img/eqlin-069.png)
![N = 7](img/eqlin-070.png)
Красным цветом показан один период оконной функции при расчете ДПФ, синим — периодические повторения, которые возникают при дискретизации ДВПФ оконной функции.
Можно видеть, что в случае периодического окна длины отсчет с индексом
не равен отсчету с индексом
. В результате при периодическом повторении мы имеем только один нулевой отсчет при стыковке, соответствующий индексу
.
Симметричное окно при имеет два крайних нулевых отсчета при
и
. В результате симметричное окно оказывается на один временно́й отсчет у́же периодического окна (можно видеть по пунктирной линии огибающей треугольного окна), что приводит к расширению основного лепестка спектральной плотности энергии окна
, как это показано на рисунке 7 для
,
и
.
![Спектральные плотности энергии периодического (красный) и симметричного (синий) окна Бартлетта для , и](img/win_periodic_freq.png)
![N = 7](img/eqlin-070.png)
![N = 8](img/eqlin-069.png)
![N = 32](img/eqlin-073.png)
Можно видеть, что для всех длин ширина основного лепестка симметричного окна больше чем у периодического. Причем можно заметить, что разница уменьшается с ростом длины
. Так для
разница ширины главного лепестка и первых боковых лепестков
становится несущественной. Однако для малых длин
потери на обработку при использовании симметричного окна будут существенно выше.
Важное свойство, которое мы можем использовать для быстрого пересчета симметричных и периодических окон можно заметить из рисунка 6. Так симметричное окно при полностью повторяет периодическое окно при
(ввиду периодичности отсчет с индексом
повторяет отсчет с индексом
). Это правило позволяет нам использовать идентичный алгоритм для расчета симметричных и периодических окон. Симметричное окно длины
получается как периодическое окно длины
плюс один симметричный отсчет равный индексу
.
В данном разделе мы приводим выражения и графики для периодических окон длины отсчетов. Симметричные окна длины
могут быть получены как периодическое окно длины
с добавлением дополнительного отсчета равного
.
Прямоугольное окно имеет единичные значения для всех :
![equation 17](img/eq-17.png)
![Прямоугольное окно](img/win_ex_rect.png)
Бартлетт использовал треугольное окно в задачах спектрального анализа:
![equation 18](img/eq-18.png)
![Треугольное окно Бартлетта](img/win_ex_bartlett.png)
Косинус окно, или синус окно, представляет собой полпериода синуса:
![equation 19](img/eq-19.png)
![Косинус окно](img/win_ex_cos.png)
Окно Ханна представляет один период приподнятого косинуса, смещенного на полдлины окна:
![equation 20](img/eq-20.png)
![Окно Ханна](img/win_ex_hann.png)
Объединив треугольное окно и приподнятый косинус получим окно Бартлетта–Ханна вида:
![equation 21](img/eq-21.png)
![Треугольное окно Бартлетта–Ханна](img/win_ex_bartlett_hann.png)
Ланцош преложил использовать главный лепесток -функции:
![equation 22](img/eq-22.png)
![Окно Ланцоша](img/win_ex_lanczos.png)
Окно Ханна использует один период приподнятого косинуса. При этом оба коэффициента (амплитуда косинуса и смещение) равны . Хемминг задал эти коэффициенты в виде:
![equation 23](img/eq-23.png)
![a_0](img/eqlin-079.png)
![equation 24](img/eq-24.png)
![Окно {Хемминга}](img/win_ex_hamming.png)
Окно Хемминга обладает хорошими частотными свойствами (низкий уровень боковых лепестков и малая ширина окна) и часто используется на практике.
Дальнейшие улучшения окон приподнятого косинуса заключалось в добавлении еще одного слагаемого. Блэкман предложил следующее окно, которое содержало три слагаемых:
![equation 25](img/eq-25.png)
![\mathrm{w}(n) = 0](img/eqlin-080.png)
![n = 0](img/eqlin-072.png)
![n = N](img/eqlin-081.png)
![\mathrm{w}(n) = 1](img/eqlin-032.png)
![n = N/2](img/eqlin-082.png)
![\alpha](img/eqlin-083.png)
![equation 26](img/eq-26.png)
![Окно Блэкмана](img/win_ex_blackman.png)
Идея добавления косинусоидальных слагаемых и оптимизации коэффициентов выглядела перспективной, в результате появилось семейство окон, содержащих уже четыре слагаемых. Одним из таких является окно Блэкмана–Харриса вида:
![equation 27](img/eq-27.png)
![equation 28](img/eq-28.png)
![Окно Блэкмана–Харриса](img/win_ex_blackman_harris.png)
Можно заметить, что уровень боковых лепестков спектральной плотности энергии уменьшается, но добавление слагаемых приводит к увеличению ширины главного лепестка.
Альтернативным окном, содержащим 4 слагаемых является окно Натталла:
![equation 27](img/eq-27.png)
![equation 29](img/eq-29.png)
![Окно Натталла](img/win_ex_nuttall.png)
Можно обратить внимание, что вид окна Натталла и его спектральная плотность энергии очень похожа на окно Блэкмана–Харриса.
Также можно отметить окно Блэкмана–Натталла, которое содержит 4 слагаемых
![equation 27](img/eq-27.png)
![equation 30](img/eq-30.png)
![Окно Блэкмана–Натталла](img/win_ex_blackman_nuttall.png)
Окно Блэкмана-Натталла обладает более низким уровнем боковых лепестков спектральной плотности мощности по сравнению с другими окнами, состоящие из четырех слагаемых косинусов.
Окно состоящие из пяти слагаемых вида:
![equation 31](img/eq-31.png)
![equation 32](img/eq-32.png)
![Окно с плоской вершиной](img/win_ex_flattop.png)
Отличие этого окна в том, что оно не ограничено значениями 0 и 1 во временно́й области, что позволяет получить широкий главный лепесток спектральной плотности энергии окна.
Данное окно использует полиномы Чебышева для представления спектра оконной функции.
Данное окно является параметрическим, и в качестве параметра окна выступает параметр максимального уровня боковых лепестков окна в частотной области, который задается в единицах децибел.
ДВПФ окна Дольф–Чебышева имеет вид:
![equation 33](img/eq-33.png)
![equation 34](img/eq-34.png)
![-\gamma](img/eqlin-085.png)
![G_C](img/eqlin-037.png)
![\gamma = 60](img/eqlin-086.png)
![-60](img/eqlin-087.png)
![G_C](img/eqlin-037.png)
Постоянный уровень боковых лепестков возможен только в случае симметричного окна. Таким образом невозможно найти периодическое окно Дольф–Чебышева. Все окна которые мы рассмотрим в данном подпараграфе являются симметричными.
Выражение для симметричного окна Дольф–Чебышева ,
имеет вид:
![equation 35](img/eq-35.png)
![equation 36](img/eq-36.png)
![equation 37](img/eq-37.png)
![equation 38](img/eq-38.png)
![\mathrm{w}(n)](img/eqlin-013.png)
![n = 0\ldots m](img/eqlin-088.png)
![equation 39](img/eq-39.png)
![Окно Дольф–Чебышева,](img/win_ex_cheby50.png)
![\gamma = 50](img/eqlin-089.png)
![Окно Дольф–Чебышева,](img/win_ex_cheby80.png)
![\gamma = 80](img/eqlin-090.png)
![Окно Дольф–Чебышева,](img/win_ex_cheby100.png)
![\gamma = 100](img/eqlin-091.png)
Использование в качестве окна колоколообразной функции приводит к семейству параметрических окон Гаусса:
![equation 40](img/eq-40.png)
Параметр позволяет менять ширину окна и получать различный уровень боковых лепестков спектральной плотности энергии оконной функции.
На рисунках показан вид окон Гаусса во временно́й и частотной областях при различном значении параметра ,
и
.
![Окно Гаусса,](img/win_ex_gaussian0p5.png)
![\alpha = 0.5](img/eqlin-092.png)
![Окно Гаусса,](img/win_ex_gaussian0p3.png)
![\alpha = 0.3](img/eqlin-093.png)
![Окно Гаусса,](img/win_ex_gaussian0p22.png)
![\alpha = 0.22](img/eqlin-094.png)
Можно заметить, что с уменьшение имеет место сужение окна во времени, расширение основного лепестка спектральной плотности энергии и одновременное уменьшение боковых лепестков.
Окно Кайзера одно из распространенных окон для синтеза цифровых КИХ-фильтров, поскольку является аппроксимацией оптимальных с точки зрения концентрации энергии в главном лепестке спектра. В этом смысле окна Кайзера позволяют получить, близкую к минимально-возможной, переходную полосу цифрового фильтра. Кроме того, возможность регулирования уровня боковых лепестков за счет изменения параметра окна, делает его весьма удобным в задачах проектирования КИХ-фильтром оконным методом.
Уравнение параметрического окна Кайзера имеет вид:
![equation 41](img/eq-41.png)
![I_0(z)](img/eqlin-095.png)
![\alpha](img/eqlin-083.png)
На рисунках 26–28 приведен вид оконной функции и ее нормированная спектральная плотность энергии.
![Окно Кайзера,](img/win_ex_kaiser4.png)
![\alpha = 4](img/eqlin-096.png)
![Окно Кайзера,](img/win_ex_kaiser8.png)
![\alpha = 8](img/eqlin-097.png)
![Окно Кайзера,](img/win_ex_kaiser12.png)
![\alpha = 12](img/eqlin-098.png)
Можно видеть, что при увеличении параметра уровень боковых лепестков оконной функции уменьшается. В практических приложениях
изменяется в пределах
.
В таблице 1 приведены основные параметры оконных функций:
- уровень максимального бокового лепестка спектральной плотности энергии
в дБ;
- нормированное когерентное усиление по амплитуде
;
- ширина главного лепестка
спектральной плотности энергии по уровню
дБ в единицах ДПФ бин;
- ширина главного лепестка
спектральной плотности энергии по уровню
дБ в единицах ДПФ бин;
- ширина главного лепестка
спектральной плотности энергии по нулевому уровню в единицах ДПФ бин;
- эквивалентная шумовая полоса
в единицах ДПФ бин;
- потери на обработку
вызванные расширением
, т.е. ухудшение отношения сигнал–шум при анализе на фоне шумов.
- максимальный уровень гребешковых искажений
дБ;
- суммарные потери на обработку в наихудшем случае, которые складываются из максимального уровня гребешковых искажений и ухудшения отношения сигнал–шум.
![Сводная таблица параметров оконных функций](img/win_table.png)
Из таблицы можно видеть, что при расширении основного лепестка спектральной плотности энергии оконной функции (увеличение ), наблюдается уменьшение гребешковых искажений
, уменьшение уровня боковых лепестков
, но вместе с тем, увеличиваются потери на обработку
ввиду увеличения
.
Таким образом, улучшая одни параметры, например увеличивая динамический диапазон спектрального анализа за счет уменьшения уровня боковых лепестков , мы ухудшаем разрешение спектрального анализа по частоте (пики гармоник становятся шире).
Кроме того при наличии шумов, мы ухудшаем отношение сигнал–шум.
Поэтому невозможно рекомендовать «универсальную оконную функцию на все случаи жизни», а нужно каждый раз, руководствуясь необходимой целесообразностью, выбирать наиболее подходящую оконную функцию, как компромисс между динамическим диапазоном, частотным разрешением, уменьшением гребешковых искажений и ухудшением отношения сигнал–шум на выходе.
В библиотеке DSPL–2.0 реализована функция window , представляющая собой коллекция оконных функций, которые были использованы для построения характеристик окон данного раздела.
Исходный код программы win_examples.c реализующий расчет характеристик рассмотренных оконных функций приведен в следующем листинге:
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
#define N 50
#define K 1000
void win_norm(double* w, int n)
{
int k;
double s = 0.0;
for(k = 0; k < n; k++)
s+=w[k];
for(k = 0; k < n; k++)
w[k] /= s;
}
void window_plot(double* w, int n, int k, int argc, char* argv[],
double ymin, double ymax,
char* png_fn, char* time_fn, char* freq_fn, char* title)
{
void* hplot; /* GNUPLOT handles */
char str[128] = {0};
double* t = NULL;
double* W = NULL;
double* f = NULL;
double fs = (double)n;
fft_t pfft = {0};
t = (double*)malloc(n* sizeof(double));
W = (double*)malloc(k*sizeof(double));
f = (double*)malloc(k*sizeof(double));
linspace(0, n, n, DSPL_PERIODIC, t);
writetxt(t, w, n, time_fn);
win_norm(w, n);
fft_mag(w, k, &pfft, fs, DSPL_FLAG_LOGMAG | DSPL_FLAG_FFT_SHIFT, W, f);
writetxt(f, W, k, freq_fn);
/* plotting 3d surface by GNUPLOT */
/* Create window 0 */
gnuplot_create(argc, argv, 820, 360, png_fn, &hplot);
gnuplot_cmd(hplot, "set grid");
memset(str,0,128);
sprintf(str, "set title '%s'", title);
gnuplot_cmd(hplot, str);
gnuplot_cmd(hplot, "set multiplot layout 1,2 rowsfirst");
gnuplot_cmd(hplot, "set xlabel 'n'");
gnuplot_cmd(hplot, "set ylabel 'w(n)'");
gnuplot_cmd(hplot, "unset key");
gnuplot_cmd(hplot, "set xrange[-4:54]");
memset(str,0,128);
sprintf(str, "set yrange[%f:%f]", ymin, ymax);
gnuplot_cmd(hplot, str);
memset(str,0,128);
sprintf(str, "plot '%s' w i lc 1,'%s' w p pt 7 ps 0.5 lc 1",
time_fn, time_fn);
gnuplot_cmd(hplot, str);
gnuplot_cmd(hplot, "set xrange[-20:20]");
gnuplot_cmd(hplot, "set yrange[-130:5]");
memset(str,0,128);
gnuplot_cmd(hplot, "set xlabel 'freq [DFT bins]'");
gnuplot_cmd(hplot, "set ylabel 'W(freq), dB'");
sprintf(str, "plot '%s' w l lc 2 ", freq_fn);
gnuplot_cmd(hplot, str);
gnuplot_close(hplot);
if(t)
free(t);
if(W)
free(W);
if(f)
free(f);
}
int main(int argc, char* argv[])
{
void* hdspl; /* DSPL handle */
hdspl = dspl_load(); /* Load DSPL function */
double w[K] = {0};
/* Rectangular window */
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_RECT, 0.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_rect.png",
"dat/win_rect_time.txt",
"dat/win_rect_freq.txt",
"Rectangular window");
/* Bartlett window (triangular)*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_BARTLETT, 0.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_bartlett.png",
"dat/win_bartlett_time.txt",
"dat/win_bartlett_freq.txt",
"Bartlett window");
/* Flat top window */
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_FLAT_TOP, 0.0);
window_plot(w, N, K, argc, argv, -1.0, 5.0,
"img/win_flattop.png",
"dat/win_flattop_time.txt",
"dat/win_flattop_freq.txt",
"Flat top window");
/* Bartlett - Hann window*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_BARTLETT_HANN, 0.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_bartletthann.png",
"dat/win_bartletthann_time.txt",
"dat/win_bartletthann_freq.txt",
"Bartlett-Hann window");
/* Blackman window*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_BLACKMAN, 0.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_blackman.png",
"dat/win_blackman_time.txt",
"dat/win_blackman_freq.txt",
"Blackman window");
/* Blackman - Harris window*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_BLACKMAN_HARRIS, 0.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_blackmanharris.png",
"dat/win_blackmanharris_time.txt",
"dat/win_blackmanharris_freq.txt",
"Blackman-Harris window");
/* Blackman - Nuttall window*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_BLACKMAN_NUTTALL, 0.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_blackmannuttall.png",
"dat/win_blackmannuttall_time.txt",
"dat/win_blackmannuttall_freq.txt",
"Blackman-Nuttull window");
/* Dolph-Chebyshev window (sidelobes level is -50 dB)*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_CHEBY, 50.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_cheby50.png",
"dat/win_cheby50_time.txt",
"dat/win_cheby50_freq.txt",
"Dolph-Chebyshev window (Rs = 50dB)");
/* Dolph-Chebyshev window (sidelobes level is -80 dB)*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_CHEBY, 80.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_cheby80.png",
"dat/win_cheby80_time.txt",
"dat/win_cheby80_freq.txt",
"Dolph-Chebyshev window (Rs = 80dB)");
/* Dolph-Chebyshev window (sidelobes level is -100 dB)*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_CHEBY, 100.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_cheby100.png",
"dat/win_cheby100_time.txt",
"dat/win_cheby100_freq.txt",
"Dolph-Chebyshev window (Rs = 100dB)");
/* Gaussian window (sigma = 0.5)*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_GAUSSIAN, 0.5);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_gaussian0p5.png",
"dat/win_gaussian0p5_time.txt",
"dat/win_gaussian0p5_freq.txt",
"Gaussian window (sigma = 0.5)");
/* Gaussian window (sigma = 0.3)*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_GAUSSIAN, 0.3);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_gaussian0p3.png",
"dat/win_gaussian0p3_time.txt",
"dat/win_gaussian0p3_freq.txt",
"Gaussian window (sigma = 0.3)");
/* Gaussian window (sigma = 0.7)*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_GAUSSIAN, 0.22);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_gaussian0p22.png",
"dat/win_gaussian0p22_time.txt",
"dat/win_gaussian0p22_freq.txt",
"Gaussian window (sigma = 0.22)");
/* Hamming window*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_HAMMING, 0.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_hamming.png",
"dat/win_hamming_time.txt",
"dat/win_hamming_freq.txt",
"Hamming window");
/* Hann window*/
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_HANN, 0.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_hann.png",
"dat/win_hann_time.txt",
"dat/win_hann_freq.txt",
"Hann window");
/* Lanczos window */
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_LANCZOS, 0.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_lanczos.png",
"dat/win_lanczos_time.txt",
"dat/win_lanczos_freq.txt",
"Lanczos window");
/* Nuttall window */
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_NUTTALL, 0.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_nuttall.png",
"dat/win_nuttall_time.txt",
"dat/win_nuttall_freq.txt",
"Nuttall window");
/* Cosine window */
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_COS, 0.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_cos.png",
"dat/win_cos_time.txt",
"dat/win_cos_freq.txt",
"Cosine window");
/* Kaiser window pi * alpha = 4 */
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_KAISER, 4.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_kaiser4p0.png",
"dat/win_kaiser4p0_time.txt",
"dat/win_kaiser4p0_freq.txt",
"Kaiser window (pi * alpha = 4)");
/* Kaiser window pi * alpha = 8 */
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_KAISER, 8.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_kaiser8p0.png",
"dat/win_kaiser8p0_time.txt",
"dat/win_kaiser8p0_freq.txt",
"Kaiser window (pi * alpha = 8)");
/* Kaiser window pi * alpha = 12 */
window(w, N, DSPL_WIN_PERIODIC | DSPL_WIN_KAISER, 12.0);
window_plot(w, N, K, argc, argv, 0.0, 1.1,
"img/win_kaiser12p0.png",
"dat/win_kaiser12p0_time.txt",
"dat/win_kaiser12p0_freq.txt",
"Kaiser window (pi * alpha = 12)");
dspl_free(hdspl); /* free dspl handle */
return 0;
}
Дискретно-временное преобразование Фурье
Дискретное преобразование Фурье
Связь ДПФ и ДВПФ. Ядро Дирихле
Индексация и перестановка спектральных отсчетов дискретного преобразования Фурье
Спектральный анализ ограниченных во времени сигналов. Эффект растекания спектра
[1] Некоторые специальные типы оконных функций обладают бо́льшими суммарными потерями, например окно с плоской вершиной, но эти окна имеют свое специальное назначение. Так окно с плоской вершиной имеет уровень гребешковых искажений всего дБ, что позволяет его использовать для прецизионного спектрального анализа, но только при высоком отношении сигнал–шум, потери которого составляют
дБ.