Спектральный анализ сигналов

image

Не так давно товарищ Makeman описывал, как с помощью спектрального анализа можно разложить некоторый звуковой сигнал на слагающие его ноты. Давайте немного абстрагируемся от звука и положим, что у нас есть некоторый оцифрованный сигнал, спектральный состав которого мы хотим определить, и достаточно точно.

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

Итак, что имеем.
Файл с отсчетами оцифрованного сигнала. Известно, что сигнал представляет собой сумму синусоид со своими частотами, амплитудами и начальными фазами, и, возможно, белый шум.

Что будем делать.
Использовать спектральный анализ для того, чтобы определить:
  • количество гармоник в составе сигнала, а для каждой: амплитуду, частоту (далее в контексте числа длин волн на длину сигнала), начальную фазу;
  • наличие/отсутствие белого шума, а при наличии, его СКО (среднеквадратическое отклонение);
  • наличие/отсутствие постоянной составляющей сигнала;
  • всё это оформить в красивенький PDF отчёт с блэкджеком и иллюстрациями.


Будем решать данную задачу на Java.

Матчасть


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

Немного особой, Фурье-магии

Не так давно, в 19 веке, французский математик Жан Батист Жозеф Фурье показал, что любую функцию, удовлетворяющую некоторым условиям (непрерывность во времени, периодичность, удовлетворение условиям Дирихле) можно разложить в ряд, который в дальнейшем получил его имя — ряд Фурье.

В инженерной практике разложение периодических функций в ряд Фурье широко используется, например, в задачах теории цепей: несинусоидальное входное воздействие раскладывают на сумму синусоидальных и рассчитывают необходимые параметры цепей, например, по методу наложения.

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

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

image
Рис. 1. Соответствие ряда Фурье и преобразования Фурье на примере амплитудного спектра.

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

Однако, преобразование Фурье сопоставляет непрерывной во времени, бесконечной функции другую, непрерывную по частоте, бесконечную функцию — спектр. Как быть, если у нас нет бесконечной во времени функции, а есть лишь какая-то записанная её дискретная во времени часть? Ответ на этот вопрос даёт дальнейшей развитие преобразования Фурье — дискретное преобразование Фурье (ДПФ).

Дискретное преобразование Фурье призвано решить проблему необходимости непрерывности и бесконечности во времени сигнала. По сути, мы полагаем, что вырезали какую-то часть бесконечного сигнала, а всю остальную временную область считаем этот сигнал нулевым.

Математически это означает, что, имея исследуемую бесконечную во времени функцию f(t), мы умножаем ее на некоторую оконную функцию w(t), которая обращается в ноль везде, кроме интересующего нас интервала времени.

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

Остальные свойства преобразования Фурье не изменяются: о них можно прочитать в соответствующей литературе.

Нам же нужно лишь знать о Фурье-образе синусоидального сигнала, который мы и будем стараться отыскать в нашем спектре. В общем случае, это пара дельта-функций, симметричная относительно нулевой частоты в частотной области.

image
Рис. 2. Амплитудный спектр синусоидального сигнала.

Я уже упомянул, что, вообще говоря, мы рассматриваем не исходную функцию, а некоторое её произведение с оконной функцией. Тогда, если спектр исходной функции — F(w), а оконной W(w), то спектром произведения будет такая неприятная операция, как свёртка этих двух спектров (F*W)(w) (Теорема о свёртке).

На практике это означает, что вместо дельта-функции, в спектре мы увидим что-то вроде этого:

image
Рис. 3. Эффект растекания спектра.

Этот эффект именуют также растеканием спектра (англ. spectral leekage). А шумы, появляющиеся вследствие растекания спектра, соответственно, боковыми лепестками (англ. sidelobes).
Для борьбы с боковыми лепестками применяют другие, непрямоугольные оконные функции. Основной характеристикой «эффективности» оконной функции является уровень боковых лепестков (дБ). Сводная таблица уровней боковых лепестков для некоторых часто используемых оконных функций приведена ниже.

Оконная функция Уровень боковых лепестков (дБ)
Окно Дирихле (прямоугольное окно) -13 дБ
Окно Ханна -31.5 дБ
Окно Хэмминга -42 дБ


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

image
Рис. 4. Отдельные спектры гармоник.

Видно, что при сложении приведённых спектров, более слабые гармоники как бы растворятся в более сильной.

image
Рис. 5. Чётко видна лишь одна гармоника. Нехорошо.

Другой подход к борьбе с растеканием спектра состоит в вычитании из сигнала гармоник, создающих это самое растекание.
То есть, установив амплитуду, частоту и начальную фазу гармоники, можно вычесть её из сигнала, при этом мы уберём и «дельта-функцию», соответствующую ей, а вместе с ней и боковые лепестки, порождаемые ей. Другой вопрос состоит в том, как же точно узнать параметры нужной гармоники. Недостаточно просто взять нужные данные из комплексной амплитуды. Комплексные амплитуды спектра сформированы по целым частотам, однако, ничто не мешает гармонике иметь и дробную частоту. В этом случае, комплексная амплитуда как бы расплывается между двумя соседними частотами, и точную её частоту, как и другие параметры, установить нельзя.

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

Посмотрим, что получится, если умножить входной сигнал на комплексную гармонику Exp(I*w*t). Спектр сигнала сдвинется на величину w вправо.
Этим свойством мы и воспользуемся, сдвигая спектр нашего сигнала вправо, до тех пор, пока гармоника не станет ещё больше напоминать дельта-функцию (то есть, пока некоторое локальное отношение сигнал/шум не достигнет максимума). Тогда мы и сможем вычислить точную частоту нужной гармоники, как w0 – wгет, и вычесть её из исходного сигнала для подавления эффекта растекания спектра.
Иллюстрация изменения спектра в зависимости от частоты гетеродина показана ниже.

image
Рис. 6. Вид амплитудного спектра в зависимости от частоты гетеродина.

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

Затем, надо оценить СКО белого шума. Хитростей здесь нет: можно просто воспользоваться формулой для вычисления СКО:

image

Автоматизируй это


Пришло время для автоматизации выделения гармоник. Повторим ещё разочек алгоритм:

1. Ищем глобальный пик амплитудного спектра, выше некоторого порога k.
1.1 Если не нашли, заканчиваем
2. Варируя частоту гетеродина, ищем такое значение частоты, при которой будет достигаться максимум некоторого локального отношения сигнал/шум в некоторой окрестности пика
3. При необходимости, округляем значения амплитуды и фазы.
4. Вычитаем из сигнала гармонику с найденной частотой, амплитудой и фазой за вычетом частоты гетеродина.
5. Переходим к пункту 1.

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

Построим функцию распределения (привет, мат. cтатистика), где по оси абсцисс будет амплитуда гармоник, а по оси ординат — количество гармоник, не превышающих по амплитуде это самое значение аргумента. Пример такой построенной функции:

image
Рис. 7. Функция распределения гармоник.

Теперь построим еще и функцию — плотность распределения. Т.е., значения конечных разностей от функции распределения.

image
Рис. 8. Плотность функции распределения гармоник.

Абсцисса максимума плотности распределения и является амплитудой гармоники, встречающейся в спектре наибольшее число раз. Отойдем от пика вправо на некоторое расстояние, и будем считать абсциссу этой точки оценкой уровня шума в нашем спектре. Вот теперь можно и автоматизировать.

Посмотреть на кусок кода, детектирующий гармоники в составе сигнала
public ArrayList<SynthesizableSignal> detectHarmonics() {
    SignalCutter cutter = new SignalCutter(source, new Signal(source));
    SynthesizableComplexExponent heterodinParameter = new SynthesizableComplexExponent();
    heterodinParameter.setProperty("frequency", 0.0);
    Signal heterodin = new Signal(source.getLength());
    Signal heterodinedSignal = new Signal(cutter.getCurrentSignal());
    Spectrum spectrum = new Spectrum(heterodinedSignal);
    int harmonic;
    while ((harmonic = spectrum.detectStrongPeak(min)) != -1) {
        if (cutter.getCuttersCount() > 10)
            throw new RuntimeException("Unable to analyze signal! Try another parameters.");
        double heterodinSelected = 0.0;
        double signalToNoise = spectrum.getRealAmplitude(harmonic) / spectrum.getAverageAmplitudeIn(harmonic, windowSize);
        for (double heterodinFrequency = -0.5; heterodinFrequency < (0.5 + heterodinAccuracy); heterodinFrequency += heterodinAccuracy) {
            heterodinParameter.setProperty("frequency", heterodinFrequency);
            heterodinParameter.synthesizeIn(heterodin);
            heterodinedSignal.set(cutter.getCurrentSignal()).multiply(heterodin);
            spectrum.recalc();
            double newSignalToNoise = spectrum.getRealAmplitude(harmonic) / spectrum.getAverageAmplitudeIn(harmonic, windowSize);
            if (newSignalToNoise > signalToNoise) {
                signalToNoise = newSignalToNoise;
                heterodinSelected = heterodinFrequency;
            }
        }
        SynthesizableCosine parameter = new SynthesizableCosine();
        heterodinParameter.setProperty("frequency", heterodinSelected);
        heterodinParameter.synthesizeIn(heterodin);
        heterodinedSignal.set(cutter.getCurrentSignal()).multiply(heterodin);
        spectrum.recalc();
        parameter.setProperty("amplitude", MathHelper.adaptiveRound(spectrum.getRealAmplitude(harmonic)));
        parameter.setProperty("frequency", harmonic - heterodinSelected);
        parameter.setProperty("phase", MathHelper.round(spectrum.getPhase(harmonic), 1));
        cutter.addSignal(parameter);
        cutter.cutNext();
        heterodinedSignal.set(cutter.getCurrentSignal());
        spectrum.recalc();
    }
    return cutter.getSignalsParameters();
}


Практическая часть


Я не претендую на звание эксперта Java, и представленное решение может быть сомнительным как по части производительности и потреблению памяти, так и в целом философии Java и философии ООП, как бы я ни старался сделать его лучше. Написано было за пару вечеров, как proof of concept. Желающие могут ознакомиться с исходным кодом на GitHub.

Единственной сложностью стала генерация PDF отчёта по результатам анализа: PDFbox ну никак не хотел работать с кириллицей. К слову, не хочет и сейчас.

В проекте использовались библиотеки:
JFreeChart – отображение графиков
PDFBox – построение отчёта
JLatexMath – рендер Latex формул

В итоге, получилась довольно массивная программа (13.6 мегабайт), удобно реализующая поставленную задачу.

Есть возможность как вырезать гармоники вручную, так и доверить эту задачу алгоритму.

Всем спасибо за внимание!


Пример отчёта, создаваемого программой.

Литература


Сергиенко А. Б. — Цифровая обработка сигналов
Поделиться публикацией
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама

Комментарии 22

    +2
    Любопытный метод, спасибо!

    А может ли он разрешить две близко лежащие синусоиды в сигнале? Допустим, мы нашли в спектре ДПФ пик (который на самом деле является «слившимися» двумя пиками от близкорасположенных гармоник). Что в этом случае даст гетеродинирование? Позволит ли оно достаточно точно определить амплитуду и фазу первой синусоиды, чтобы, после вычитания ее из сигнала, алгоритм смог детектировать вторую синусоиду?
      +1
      Скорее нет, чем да. Насколько я понял, речь идёт о близости гармоник в диапазоне ± 1 единицы частоты (опять же, в контексте числа длин волн на сигнал). Гетеродинирование не решает проблему недостатка разрешающей способности по частоте. Увы. Правда, при определенной разнице между частотами гармоник, можно, по крайней мере приблизительно, определить соотношение амплитуд, и, собственно, то, что гармоник — две.
        +1
        Я вот читал в книге С.Л. Марпл-Мл. «Цифровой спектральный анализ и его приложения» о модифицированном методе Прони, который пытается представить фрагмент сигнала в виде суммы незатухающих синусоид произвольной частоты. Соответственно, спектр по мод. Прони является строго линейчатым. Там даже были какие-то шансы разделить близкорасположенные гармоники лучше, чем это возможно с помощью ДПФ. Правда, я этот метод так и не реализовал и не испытал.

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

        Также рекомендую более подходящую, чем Java, платформу для отладки и испытания методов обработки сигналов — Matlab.
          0
          Очень удобно для испытания методов обработки сигналов использовать Python с библиотекой SciPy как бесплатная альтернатива Matlab.
            0
            Любопытно. Спасибо, отправил в список на прочтение.

            Сам пользуюсь лицензионной Wolfram Mathematica для подтверждения принципа. Дело привычки наверно. В данном случае Java выбрал можно сказать случайно. Устал от веб-разработки, решил что-нибудь standalone написать.
        0
        Фраза «Но как показал наш спектральный анализ...» больше не вызывает у меня улыбку.
          0
          Вопрос к формуле с рассчетом СКО, если мне не изменяет память, в таком случае она получится смещенной.
            0
            Да, не изменяет. Оценка СКО будет смещённой. Намекаете на поправочный множитель n/(n-1)?
              +1
              Да вроде прямым текстом говорю. Не знаю, после универа деление для дисперсии на n-1 в подкорку въедается уже, прям глаз режет, когда не так. При больших n конечно разница будет не такой уж большой, но все-таки математически это неверно.
                0
                Принимается. Спасибо за замечание! В ближайшее время поправлю.
                  0
                  Вот что говорят по этому поводу William H. Press et al в книге «Numerical recipes in C++»:

                  «We might also comment that if the difference between N and N-1 ever matters to you, then you are probably up to no good anyway — e.g., trying to substantiate a questionable hypothesis with marginal data».
                    0
                    Чего еще ожидать от физика-практика. Если оценка смещенная, в общем случае ей нельзя пользоваться. И если он утверждает «а, пофиг, все равно разница небольшая», то это говорит не в его пользу. У меня в универе тоже есть на кафедре куча заслуженных доцентов, которые всякий маразматичный бред несут, хотя в свое время были очень крутыми ребятами/девчатами. Если человек застрял в изучении случайных процессов на уровне учебника Вентцеля — то для студента 2-3 курса техвуза это очень даже неплохо, но вот для подобных утверждений явно недостаточно.

                    Ничего не хочу плохого сказать про сего достопочтенного мэтра, но в данном случае он не прав. Очень часто N не такое уж большое число, а при малых N вклад будет очень существенным. А заниженная оценка дисперсии это неприятная штука.

                    Я сейчас пишу черновик на тему регрессионного/дисперсионного анализов, как выложу, сможете оценить разницу при таком «наплевательском» отношении к константам. Вычесть единичку это очень трудно, я понимаю, но всё же. Ну и нужно понимать, что если такая небрежность наблюдается по всей программе, то ошибка будет накапливаться, что в результате даст ощутимую погрешность, а то и решающую.
                      +1
                      Вот было бы неплохо привести пример, при каких условиях пренебрежение N-1 может дать ощутимую или решающую погрешность в результате обработки данных. Потому что математика — это одно, а физика вносит свои коррективы.

                      А то при чисто «математическом» подходе, если вы измерили массу три раза и получили 100, 100 и 101кг — вы можете посчитать среднее, которое равно 100,333333333333кг с точностью ~1e-15. С такой точностью не известна даже масса национальных эталонов килограмма. Но повода для гордости здесь нет. Какой смысл считать 15 цифр после запятой, если погрешность измерений составляет порядка 1e-2?

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

                        Не вопрос: допустим я хочу понять зависимость от каких-нибудь факторов по странам мира. Это всего лишь N=200, соответственно относительная погрешность будет 0.5% при подсчете одной только дисперсии, при использовании в дальнейшем ошибка будет накапливаться, т.к. не в одной только дисперсии забьют на правильную формулу. Если же буду брать статистику по странам Европы, то получу всего несколько десятков точек. При активном использовании в формулах ошибка может достигнуть 30-40% на ровном месте. Из-за того, что кому-то лень написать единичку.

                        Поделили на n — получили заниженную оценку дисперсии. На практике дисперсия будет больше. Скорее всего ненамного, но это всё относительно. И при малых n может быть крайне существенно.
              +1
              Несколько месяцев назад мы с коллегой столкнулись именно с рассмотренной задачей, только формулировка была несколько иная: есть измерительное устройство, которое регистрирует и оцифровывает синусоидальный сигнал с шумами, и в микроконтроллере нужно было точно рассчитать, какое кол-во оцифрованных отсчетов приходится на один период преобладающей гармоники. Понятно, что в общем случае кол-во отсчетов в периоде должно получиться дробным, а точность требовалась аж до 5-го знака после запятой. Идея алгоритма, который я применил, была схожа с Вашим решением, но только преобразование Фурье использовать не пришлось. Я просто вычислял несколько «чистых» синусоид, имеющих частоту из окрестности поиска, и считал свертку с оцифрованным сигналом — это фактически была корреляционная функция, и ее экстремум соответствовал искомому значению. Однако кучу раз считать синусы в МК дело затратное, поэтому решено было «чистые» синусы заменить на периодическую кусочно-линейную «пилу», от чего точность совсем не пострадала.
              Так что с практической точки зрения полезная статья :)
                0
                Рад, что вам понравилось ) Интересная идея. И как в плане производительности на МК оно?
                  +1
                  Для выборки сигнала, представленной более чем 3000 отсчетами, занимает около нескольких секунд (3-4). Но это, разумеется, зависит от того, какая арифметика используется — целочисленная либо с плавающей точкой (для integer производительность раза в 1.5 больше). Микроконтроллер из семейства stm32.
                  0
                  Я бы вам рекомендовал для указанных целей попробовать систему ФАПЧ (в английской терминологии — PLL). Доказано, что эта система обеспечивает максимальное соотношение «точность / время измерения» при измерении частоты. К тому же, она мало затратная в реализации. Матчасть только надо изучить, там она нетривиальная.
                  0
                  Подскажите, пожалуйста, почему по оси амплитуд на диаграммах амплитудного спектра значение этих самых амплитуд как-то, мягко говоря, не очень похоже на амплитуды, которые содержит исходный сигнал? Как правило, амплитуда на амплитудном спектре существенно (на порядки) выше амплитуды исходного сигнала, если смотреть на pdf с примерами.
                  Я полагаю, здесь достаточно простой ответ из области математики. Просто я от этой области далек. Может кто подскажет?
                  Автору спасибо за статью!
                    0
                    Это не баг, это фича. Амплитуду из спектра следует поделить на половину длины сигнала. А для нулевой частоты — на всю длину сигнала.
                      0
                      Спасибо за ответ!
                      Длина сигнала в данном случае — это в секундах или в числе дискретных отрезков?
                        0
                        Длина сигнала в данном контексте — число его (сигнала) отсчётов.

                  Только полноправные пользователи могут оставлять комментарии. Войдите, пожалуйста.

                  Самое читаемое