Функции оконного сглаживания

Содержание

DSPL-2.0 — свободная библиотека алгоритмов цифровой обработки сигналов

Распространяется под лицензией LGPL v3

Страница проекта на SourceForge libdspl-2.0 on SourceForge

Обнаружили ошибку? Выделите ее мышью и нажмите ctrl+enter
ДВПФ гармонического сигнала, ограниченного по времени оконной функцией

В предыдущем разделе мы анализировали эффект растекания спектра дискретного преобразования Фурье (ДПФ) и говорили, что он возникает в результате дискретизации дискретно-временного преобразования Фурье (ДВПФ) ограниченного во времени сигнала.

Пара ДВПФ сигнала s(n), n = 0\ldots N-1 имеет вид:

equation 1
(1)

ДВПФ спектр дискретного сигнала представляет собой 2\pi-периодическую функцию нормированной частоты \omega.

Взяв N дискретных значение S(\omega) на сетке частот \omega_k = 2\pi k/ N,  k = 0\ldots N-1 получим пару дискретного преобразования Фурье:

equation 2
(2)

Рассмотрим дискретно-временное преобразование Фурье комплексного гармонического сигнала s(n) = \exp(j \omega_0 n), ограниченной длительности N отсчетов, частота \omega_0 изменяется в пределах от -\pi до \pi. Ограничение по длительности сигнала s(n) осуществляется умножением на оконную функцию \mathrm{w}(n). Тогда ДВПФ S(\omega) равно:

equation 3
(3)

Таким образом, ДВПФ комплексного гармонического сигнала полностью описывается ДВПФ оконной функции, которая используется для ограничения длительности сигнала.

В предыдущем разделе мы рассмотрели эффект растекания ДПФ, который возникает, когда частота гармонического сигнала \omega_0 не совпадает с частотами дискретизации ДВПФ. Мы отмечали, что эффект растекания спектра приводит к гребешковым искажениям амплитуды гармонического сигнала, а также к уменьшению динамического диапазона спектрального анализа ввиду маскирования слабых сигналов под боковыми лепестками мощных гармоник при растекании спектра последних.

Оконное сглаживание как способ уменьшения эффекта растекания спектра

Эффект растекания ДПФ возникает в результате дискретизации ДВПФ ограниченного во времени сигнала на боковых лепестках спектра оконной функции W(\omega). Для того чтобы уменьшить растекание, надо уменьшить уровень боковых лепестков спектра W(\omega).

Как мы знаем, уровень боковых лепестков спектра и скорость их убывания по частоте зависит от степени гладкости исходной функции. Поэтому для уменьшения эффекта растекания необходимо изменить оконную функцию  \mathrm{w}(n), исключив скачки в начале и конце. На рисунке показаны \mathrm{w}(n) в виде прямоугольного окна (несглаженная) и сглаженное окно Хэмминга, а также их спектральные плотности энергии |W(\omega)|^2.

ДВПФ прямоугольного окна и окна Хемминга,
Рисунок 1. ДВПФ прямоугольного окна и окна Хемминга, N=50

Можно видеть, что уровень боковых лепестков спектра окна Хемминга значительно ниже чем у прямоугольного окна. Таким образом, использовав сглаженные окна, можно получить увеличение динамического диапазона спектрального анализа, и также уменьшить уровень гребешковых искажений.

Пересчет единиц частоты спектрального анализа

В спектральном анализе спектры ДВПФ и ДПФ могут выводится на график и представляется в различных единицах частоты. Например ДВПФ S(\omega) определяется как 2\pi периодическая функция нормированной частоты \omega. Тогда один период ДВПФ может быть переставлен в диапазоне от 0 дo 2\pi. В этом случае, ввиду периодичности S(\omega), частотный диапазон от \pi до 2 \pi соответствует отрицательной области от -\pi до нуля, как это показано на рисунке 2а.

Если для вывода спектра на график необходимо \omega = 0 разместить в центре, то можно представить частоту \omega от -\pi до \pi и произвести перестановку спектральных компонент, как это показано на рисунке 2б.

Пересчет единиц частоты
Рисунок 2. Пересчет единиц частоты

От нормированной частоты ДВПФ \omega можно перейти к спектральному представлению S(f), где f — частота в единицах герц. Для этого необходимо знать частоту дискретизации исходного сигнала F_{\text{s}}, которая использовалась для оцифровки данных. Тогда f = \omega F_{\text{s}}  / (2 \pi) Гц, как это показано на рисунках 2б и г.

Ну и наконец, для анализа характеристик оконной функции мы будем нормировать частоту \omega к частоте f в единицах спектральных отсчетов ДПФ (ДПФ бин). В этом случае нормировка задается выражением f = N F_{\text{s}}  / (2 \pi) ДПФ бин, где N — размер выборки данных (количество точек ДПФ). Тогда f=1 ДПФ бин соответсвует первому отсчету ДПФ S(1). Графически шкала частот в единицах ДПФ бин показана на рисунках 2д и е.

Основные параметры оконных функций

Можно многими способами сгладить прямоугольное окно, и таким образом, получить множество семейств различных оконных функций. Поэтому нам надо рассмотреть параметры оконных функций, которые так или иначе влияют на качество спектрального анализа. Имея набор параметров мы сможем сравнить окна между собой и получить объективный подход для выбора того или иного окна для практического применения.

Когерентное усиление

При совпадении частоты \omega_0 с частотой некоторого бина ДПФ имеет место соотношение \omega_0 = 2 \pi m /N, где m — цело число от 0 до N-1. Тогда ДПФ S(m) = W(0) и имеет максимально-возможное значение. В случае прямоугольного окна \mathrm{w}(n) = 1 при n = 0\ldots N-1, тогда W(0) = N и имеет место максимальное когерентное усиление равное N среди всех возможных окон. В этом случае комплексная экспонента соответствующая m-му бину ДПФ представляет собой согласованный фильтр, для гармонического сигнала на частоте \omega_0 = 2 \pi m /N и мы имеем максимальные отклик на выходе согласованного фильтра. Таким образом, величина когерентного усиления на выходе ДПФ для комплексной экспоненты соответствующей m-му бину ДПФ равно G_C = |W(0)| раз по амплитуде или

equation 4
(4)
Максимальное когерентное усиление G_C = N раз будет для прямоугольного окна. Использование сглаженного непрямоугольного окна будет уменьшать величину отклика на выходе в |W(0)|/N раз по амплитуде, или на величину

equation 5
(5)
относительно отклика с прямоугольным окном.

Еще раз заметим, что когерентное усиление G_C относится только к случаю, когда частота гармонического сигнала попадает на бин ДПФ и эффекта растекания не наблюдается.

Уровень гребешковых искажений

Если частота гармонического сигнала s(n) = \exp(j \omega_0 n) не совпадает с частотой бина ДПФ, то наблюдаются гребешковые искажения амплитуды гармонического сигнала на выходе ДПФ, ввиду того, что дискретизация ДВПФ производится на боковых лепестках спектра оконной функции W(2\pi k / N - \omega_0).

Поскольку ДВПФ оконной функции W(\omega) является симметричной функцией относительно нулевой частоты, то смещение частоты \omega_0 на половину бина ДПФ приведет к максимальному уровню гребешковых искажений L_S, соответствующему частоте  \omega  = 0.5 \cdot 2 \pi / N = \pi /N.

Для исключения зависимости от длины выборки N можно нормировать уровень гребешковых искажений к уровню когерентного усиления. Тогда

equation 6
(6)
или

equation 7
(7)

Спектральные характеристики оконных функций

На рисунке 3 показаны частотные параметры спектральной плотности энергии |W(f)|^2 оконных функций.

Спектральные характеристики оконной функции
Рисунок 3. Спектральные характеристики оконной функции

Значение частоты f на рисунке 3 приведено в единицах спектральных отсчетов ДПФ (бин ДПФ).

На графиках можно видеть значение когерентного усиления G_C (в данном случае усиление не нормировано к значению усиления прямоугольного окна). Относительно G_C можно выделить уровень максимального бокового лепестка \gamma_{SL} оконной функции.

Главный лепесток оконной функции характеризуется шириной по уровню -3 дБ и -6 дБ (полосы \textrm{BW}_{-3 \textrm{ дБ}} и \textrm{BW}_{-6 \textrm{ дБ}} соответственно), а также шириной главного лепестка по нулевому уровню \textrm{BW}_{0}.

Значение |W(f)|^2 на частоте f = \pm0.5 бина ДПФ показывает максимальный уровень гребешковых искажений L_S относительно уровня G_C.

Эквивалентная шумовая полоса
Рассмотрим квадрат модуля ДВПФ W(\omega) оконной функции \mathrm{w}(n) длительности N отсчетов:

equation 8
(8)
где \omega — нормированная круговая частота. Функция W(\omega) является 2\pi-периодической функцией, тогда обратное ДВПФ имеет вид:

equation 9
(9)
Эквивалентной шумовой полосой оконной функции называется такая полоса пропускания \mathrm{ENBW} идеального прямоугольного фильтра, которая обеспечивает одинаковую дисперсию белого гауссова шума на выходе, при равном коэффициенте передачи на нулевой частоте. В случае оконного сглаживания, \mathrm{ENBW} показывает полосу частот в которой заключена таже энергия, что и энергия оконной функции со спектральной плотностью энергии |W(\omega)|^2 при равном значении |W(0)|^2, как это показно на рисунке 4.

Эквивалентная шумовая полоса оконной функции
Рисунок 4. Эквивалентная шумовая полоса оконной функции
Таким образом эквивалентную шумовую полосу можно рассчитать из уравнения:

equation 10
(10)
Из уравнения (8):

equation 11
(11)
Рассмотрим теперь энергию окна во временно́й области:

equation 12
(12)
Подставим в (12) выражение обратного ДВПФ (9):

equation 13
(13)

Тогда можно заключить, что

equation 14
(14)
откуда с учетом (11)

equation 15
(15)
Для выражения \mathrm{ENBW} в единицах отсчетов ДПФ необходимо вспомнить, что \omega = 0 соответствует нулевому отсчету ДПФ, а \omega = 2\pi соответствует бину ДПФ с номером N. Тогда

equation 16
(16)
Заметим, что для прямоугольного окна \mathrm{w}(n) = 1, n = 0\ldots N-1, тогда \mathrm{ENBW}= 1 бин ДПФ. Для сглаженных окон, данный параметр всегда больше единицы, что объясняется более широким основным лепестком спектральной плотности энергии |W(\omega)|^2.

Потери на обработку. Ухудшение отношения сигнал–шум
С шириной эквивалентной шумовой полосы связан параметр L_P, который характеризует ухудшение отношения сигнал–шум при анализе гармонического сигнала на частоте бина ДПФ на фоне белого гауссова шума.

Поскольку эквивалентная шумовая полоса сглаженных окон выше чем прямоугольного окна, то это означает более высокий уровень шума на выходе.

Это можно пояснить так. Поскольку прямоугольное окно обладает максимальным когерентным усилением и является согласованным фильтром для гармонического сигнала на частоте бина ДПФ, то оно обеспечивает максимальное отношение сигнал–шум на выходе. Любое другое сглаженное окно будет обеспечивать худшее отношение сигнал–шум для гармонического сигнала на частоте бина ДПФ.

На рисунке 5 показан результат оценки спектральной плотности мощности сигнала s(n), содержащего две гармоники на частотах f_0 = 100 Гц и f_1 = 385 Гц на фоне белого гауссова шума при частоте дискретизации 1 кГц и использовании ДПФ размера N = 100 точек и различных оконных функциях. Оценка велась при усреднении по ансамблю из 10000 реализаций.

Потери при обработке гармонических сигналов
Рисунок 5. Потери при обработке гармонических сигналов

Заметим, что при построении графиков, разница в когерентном усилении оконных функций была скомпенсирована и выровнена с усилением прямоугольного окна. В результате значение гармоник на частоте 100 Гц для всех окон равно 20 дБ.

Можно видеть, что прямоугольное окно обеспечивает максимальное отношение сигнал–шум D = 19 дБ для сигнала на частоте 100 Гц, но при этом наблюдаются максимальные гребешковые искажения L_S = 3.94 дБ.

При использовании окна Хемминга величина гребешковых искажений уменьшилась до значения L_S = 1.75 дБ, однако отношение сигнал–шум также ухудшилось на величину L_P = 1.34 дБ.

Для окна Блэкмана-Харриса гребешковые искажения составляют L_S = 0.83 дБ, в то время как отношение сигнал–шум ухудшается на величину L_P = 3.02 дБ.

Можно обратить внимание, что суммарные потери на обработку, которые включают гребешковые искажения и ухудшения отношения сигнал–шум для наихудшего случая, т.е. сигнала на частоте 385 Гц составляют L_S = 3.94 дБ для прямоугольного окна, L_S+L_P = 3.09 дБ для окна Хемминга и L_S+L_P = 3.85 дБ для окна Блэкмана-Харриса. Таким образом, в наихудшем случае, суммарные потери на обработку при использовании прямоугольного окна, как правило, остаются выше, чем при использование сглаженных окон[1].

Данный пример был рассчитан с использованием библиотеки DSPL–2.0. Исходный код программы scalloping_loss.c расчета данных для построения графиков рисунка 5 приведён в следующем листинге.


#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;
}

Симметричные и периодические оконные функции

Оконные функции используются в различных задачах спектрального анализа, а также при синтезе цифровых КИХ-фильтров. В задачах спектрального анализа мы стремимся получить более узкий основной лепесток спектра оконной функции, в то время как в задачах синтеза КИХ-фильтров часто требуется обеспечить линейность ФЧХ, что требует симметрии оконной функции во временно́й области.

Для удовлетворения требований спектрального анализа и задач синтеза линейнофазовых КИХ-фильтров различают периодические и симметричные оконные функции. Периодические оконные функции обладают лучшими частотными свойствами, но не имеют свойства временно́й симметрии. Пример периодических и симметричных треугольных окон Бартлетта показан на рисунке 6 для четной N=8 и нечетной N=7 длин окна.

Периодические и симметричные окна Бартлетта для  и
Рисунок 6. Периодические и симметричные окна Бартлетта для N = 8 и N = 7

Красным цветом показан один период оконной функции при расчете ДПФ, синим — периодические повторения, которые возникают при дискретизации ДВПФ оконной функции.

Можно видеть, что в случае периодического окна длины N = 8 отсчет с индексом n = 7 не равен отсчету с индексом n = 0. В результате при периодическом повторении мы имеем только один нулевой отсчет при стыковке, соответствующий индексу n = 0.

Симметричное окно при N=8 имеет два крайних нулевых отсчета при n = 0 и n = 7. В результате симметричное окно оказывается на один временно́й отсчет у́же периодического окна (можно видеть по пунктирной линии огибающей треугольного окна), что приводит к расширению основного лепестка спектральной плотности энергии окна |W(f)|^2, как это показано на рисунке 7 для N = 7, N = 8 и N = 32.

Спектральные плотности энергии периодического (красный) и симметричного (синий) окна Бартлетта для ,  и
Рисунок 7. Спектральные плотности энергии периодического (красный) и симметричного (синий) окна Бартлетта для N = 7, N = 8 и N = 32

Можно видеть, что для всех длин N ширина основного лепестка симметричного окна больше чем у периодического. Причем можно заметить, что разница уменьшается с ростом длины N. Так для N = 32 разница ширины главного лепестка и первых боковых лепестков |W(\omega)|^2 становится несущественной. Однако для малых длин N потери на обработку при использовании симметричного окна будут существенно выше.

Важное свойство, которое мы можем использовать для быстрого пересчета симметричных и периодических окон можно заметить из рисунка 6. Так симметричное окно при N=8 полностью повторяет периодическое окно при N = 7 (ввиду периодичности отсчет с индексом n= 7 повторяет отсчет с индексом n =0). Это правило позволяет нам использовать идентичный алгоритм для расчета симметричных и периодических окон. Симметричное окно длины N получается как периодическое окно длины N-1 плюс один симметричный отсчет равный индексу n = 0.

Некоторые широко применяемы типы оконных функций

В данном разделе мы приводим выражения и графики для периодических окон длины N отсчетов. Симметричные окна длины N могут быть получены как периодическое окно длины N-1 с добавлением дополнительного отсчета равного \mathrm{w}(0).

Прямоугольное окно

Прямоугольное окно имеет единичные значения для всех n = 0\ldots N-1:

equation 17
(17)
Вид прямоугольного окна и его спектральная плотность энергии приведены на рисунке 8

Прямоугольное окно
Рисунок 8. Прямоугольное окно

Треугольное окно Бартлетта

Бартлетт использовал треугольное окно в задачах спектрального анализа:

equation 18
(18)
Вид треугольного окна Бартлетта и его спектральная плотность энергии приведены на рисунке 9

Треугольное окно Бартлетта
Рисунок 9. Треугольное окно Бартлетта

Косинус окно (синус окно)

Косинус окно, или синус окно, представляет собой полпериода синуса:

equation 19
(19)
Вид косинус окно и его спектральная плотность энергии приведены на рисунке 10
Косинус окно
Рисунок 10. Косинус окно

Окно Ханна

Окно Ханна представляет один период приподнятого косинуса, смещенного на полдлины окна:

equation 20
(20)
Вид треугольного окна Ханна и его спектральная плотность энергии приведены на рисунке 11
Окно Ханна
Рисунок 11. Окно Ханна

Окно Бартлетта-Ханна

Объединив треугольное окно и приподнятый косинус получим окно Бартлетта–Ханна вида:

equation 21
(21)
Вид треугольного окна Ханна и его спектральная плотность энергии приведены на рисунке 12.

Треугольное окно Бартлетта–Ханна
Рисунок 12. Треугольное окно Бартлетта–Ханна

Окно Ланцоша

Ланцош преложил использовать главный лепесток \operatorname{sinc}-функции:

equation 22
(22)
Вид окна Ланцоша и его спектральная плотность энергии приведены на рисунке 13.
Окно Ланцоша
Рисунок 13. Окно Ланцоша

Окно Хемминга

Окно Ханна использует один период приподнятого косинуса. При этом оба коэффициента (амплитуда косинуса и смещение) равны 0.5. Хемминг задал эти коэффициенты в виде:

equation 23
(23)
и оптимизировал a_0 с точки зрения частотных свойств окна. В результате получил следующую оконную функцию:

equation 24
(24)
Окно {Хемминга}
Рисунок 14. Окно {Хемминга}

Окно Хемминга обладает хорошими частотными свойствами (низкий уровень боковых лепестков и малая ширина окна) и часто используется на практике.

Окно Блэкмана

Дальнейшие улучшения окон приподнятого косинуса заключалось в добавлении еще одного слагаемого. Блэкман предложил следующее окно, которое содержало три слагаемых:

equation 25
(25)
Можно заметить, что \mathrm{w}(n) = 0 при n = 0 и n = N, а также \mathrm{w}(n) = 1 при n = N/2. Оптимизируя параметр \alpha было найдено следующее окно:

equation 26
(26)
Окно содержащее три слагаемых обладает более низким уровнем боковых лепестков спектральной плотности энергии, как это показано на рисунке 15

Окно Блэкмана
Рисунок 15. Окно Блэкмана

Окно Блэкмана-Харриса

Идея добавления косинусоидальных слагаемых и оптимизации коэффициентов выглядела перспективной, в результате появилось семейство окон, содержащих уже четыре слагаемых. Одним из таких является окно Блэкмана–Харриса вида:

equation 27
(27)
equation 28
(28)
Вид окна Блэкмана–Харриса и его спектральная плотность энергии приведены на рисунке 16.

Окно Блэкмана–Харриса
Рисунок 16. Окно Блэкмана–Харриса

Можно заметить, что уровень боковых лепестков спектральной плотности энергии уменьшается, но добавление слагаемых приводит к увеличению ширины главного лепестка.

Окно Натталла

Альтернативным окном, содержащим 4 слагаемых является окно Натталла:

equation 27
(27)
equation 29
(29)
Вид окна Натталла и его спектральная плотность энергии приведены на рисунке 17.
Окно Натталла
Рисунок 17. Окно Натталла

Можно обратить внимание, что вид окна Натталла и его спектральная плотность энергии очень похожа на окно Блэкмана–Харриса.

Окно Блэкмана-Натталла

Также можно отметить окно Блэкмана–Натталла, которое содержит 4 слагаемых

equation 27
(27)
equation 30
(30)
Вид окна Натталла и его спектральная плотность энергии приведены на рисунке 18.

Окно Блэкмана–Натталла
Рисунок 18. Окно Блэкмана–Натталла

Окно Блэкмана-Натталла обладает более низким уровнем боковых лепестков спектральной плотности мощности по сравнению с другими окнами, состоящие из четырех слагаемых косинусов.

Окно с плоской вершиной

Окно состоящие из пяти слагаемых вида:

equation 31
(31)
equation 32
(32)
обладает плоской вершиной в частотной области, как это показано на рисунке 19.

Окно с плоской вершиной
Рисунок 19. Окно с плоской вершиной

Отличие этого окна в том, что оно не ограничено значениями 0 и 1 во временно́й области, что позволяет получить широкий главный лепесток спектральной плотности энергии окна.

Параметрическое окно Дольф–Чебышева

Данное окно использует полиномы Чебышева для представления спектра оконной функции. Данное окно является параметрическим, и в качестве параметра окна выступает параметр \gamma максимального уровня боковых лепестков окна в частотной области, который задается в единицах децибел. ДВПФ окна Дольф–Чебышева имеет вид:

equation 33
(33)
equation 34
(34)
Окно Дольф–Чебышева имеет постоянный уровень боковых лепестков в частотной области на уровне -\gamma дБ, относительно уровня когерентного усиления G_C. Так параметр окна \gamma = 60 означает, что в результате синтеза будет получим окно с постоянным уровнем боковых лепестков -60 дБ относительно G_C.

Постоянный уровень боковых лепестков возможен только в случае симметричного окна. Таким образом невозможно найти периодическое окно Дольф–Чебышева. Все окна которые мы рассмотрим в данном подпараграфе являются симметричными.

Выражение для симметричного окна Дольф–Чебышева \mathrm{w}(n), n = 0\ldots N-1 имеет вид:

equation 35
(35)
equation 36
(36)
equation 37
(37)
equation 38
(38)
Для расчета симметричного окна достаточно рассчитать одну половину \mathrm{w}(n) для n = 0\ldots m, а вторую половину можно скопировать:

equation 39
(39)

Окно Дольф–Чебышева,
Рисунок 20. Окно Дольф–Чебышева, \gamma = 50

Окно Дольф–Чебышева,
Рисунок 21. Окно Дольф–Чебышева, \gamma = 80

Окно Дольф–Чебышева,
Рисунок 22. Окно Дольф–Чебышева, \gamma = 100

Параметрическое окно Гаусса

Использование в качестве окна колоколообразной функции приводит к семейству параметрических окон Гаусса:

equation 40
(40)

Параметр \alpha позволяет менять ширину окна и получать различный уровень боковых лепестков спектральной плотности энергии оконной функции.

На рисунках показан вид окон Гаусса во временно́й и частотной областях при различном значении параметра \alpha = 0.5, \alpha = 0.3 и \alpha = 0.22.

Окно Гаусса,
Рисунок 23. Окно Гаусса, \alpha = 0.5

Окно Гаусса,
Рисунок 24. Окно Гаусса, \alpha = 0.3

Окно Гаусса,
Рисунок 25. Окно Гаусса, \alpha = 0.22

Можно заметить, что с уменьшение \alpha имеет место сужение окна во времени, расширение основного лепестка спектральной плотности энергии и одновременное уменьшение боковых лепестков.

Параметрическое окно Кайзера

Окно Кайзера одно из распространенных окон для синтеза цифровых КИХ-фильтров, поскольку является аппроксимацией оптимальных с точки зрения концентрации энергии в главном лепестке спектра. В этом смысле окна Кайзера позволяют получить, близкую к минимально-возможной, переходную полосу цифрового фильтра. Кроме того, возможность регулирования уровня боковых лепестков за счет изменения параметра окна, делает его весьма удобным в задачах проектирования КИХ-фильтром оконным методом.

Уравнение параметрического окна Кайзера имеет вид:

equation 41
(41)
где I_0(z) — модифицированная функция Бесселя первого рода нулевого порядка, \alpha — параметр окна.

На рисунках 26–28 приведен вид оконной функции и ее нормированная спектральная плотность энергии.

Окно Кайзера,
Рисунок 26. Окно Кайзера, \alpha = 4

Окно Кайзера,
Рисунок 27. Окно Кайзера, \alpha = 8

Окно Кайзера,
Рисунок 28. Окно Кайзера, \alpha = 12

Можно видеть, что при увеличении параметра \alpha уровень боковых лепестков оконной функции уменьшается. В практических приложениях \alpha изменяется в пределах 2 < \alpha < 16.

Сводная таблица параметров оконных функций

В таблице 1 приведены основные параметры оконных функций:

Сводная таблица параметров оконных функций
Таблица 1. Сводная таблица параметров оконных функций

Из таблицы можно видеть, что при расширении основного лепестка спектральной плотности энергии оконной функции (увеличение \mathrm{BW}_0), наблюдается уменьшение гребешковых искажений L_S, уменьшение уровня боковых лепестков \gamma_{SL}, но вместе с тем, увеличиваются потери на обработку L_P ввиду увеличения \mathrm{ENBW}.

Таким образом, улучшая одни параметры, например увеличивая динамический диапазон спектрального анализа за счет уменьшения уровня боковых лепестков \gamma_{SL}, мы ухудшаем разрешение спектрального анализа по частоте (пики гармоник становятся шире). Кроме того при наличии шумов, мы ухудшаем отношение сигнал–шум.

Поэтому невозможно рекомендовать «универсальную оконную функцию на все случаи жизни», а нужно каждый раз, руководствуясь необходимой целесообразностью, выбирать наиболее подходящую оконную функцию, как компромисс между динамическим диапазоном, частотным разрешением, уменьшением гребешковых искажений и ухудшением отношения сигнал–шум на выходе.

Реализация в DSPL--2.0

В библиотеке 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] Некоторые специальные типы оконных функций обладают бо́льшими суммарными потерями, например окно с плоской вершиной, но эти окна имеют свое специальное назначение. Так окно с плоской вершиной имеет уровень гребешковых искажений всего L_S = 0.02 дБ, что позволяет его использовать для прецизионного спектрального анализа, но только при высоком отношении сигнал–шум, потери которого составляют L_P = 5.76 дБ.

Список литературы
[1] Harris F. On the Use of Windows for Harmonic Analysis With the Discrete Fourier Transform Proceedings of the IEEE, Jan. 1978, Vol. 66, Num. 1, 51–83

[2] Оппенгейм А., Шаффер Р. Цифровая обработка сигналов. Москва, Техносфера, 2012. 1048 с. ISBN 978-5-94836-329-5

[3] Сергиенко А.Б. Цифровая обработка сигналов СПб, Питер, 2002.

[4] Рабинер, Л., Гоулд, Б. Теория и применение цифровой обработки сигналов. Москва, Мир, 1978. 848 с.

[5] Марпл-мл., С.Л. Цифровой спектральный анализ и его приложения. Москва, Мир, 1990, 584 c.

[6] Лайонс Р. Цифровая обработка сигналов. Москва, ООО «Бином-Пресс», 2006, 656 c.

[7] Гольденберг Л.М., Матюшкин Б.Д., Поляк М.Н. Цифровая обработка сигналов. Москва, Радио и связь, 1985.

Последнее изменение страницы: 20.01.2024 (19:29:12)
Страница создана Latex to HTML translator ver. 5.20.11.14