Генератор сигналов с визуализацией спектра

[The trick]

Всем привет!

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

1. Белый шум;
2. Розовый шум;
3. Коричневый/красный шум;
4. Синий шум;
5. Фиолетовый шум;
6. Синусоида с изменяющейся частотой;
7. Прямоугольная волна (нечетные гармоники убывающие на 6db/oct);
8. Пилообразная волна (все гармоники убывающие на 6db/oct);
9. Треугольная волна (нечетные гармоники убывающие на 12db/oct);
10. Сигнал только с четными гармониками;
11. Сложная частотная модуляция;
12. Амплитудная модуляция;

Визуализатор спектра поддерживает разрешение до 32768 точек, а также усреднение спектра по нескольким выборкам. Визуализатор оформлен в виде отдельного контрола и может быть использован в других проектах. Для визуализации спектра применяется быстрое преобразование Фурье (FFT). Данная процедура оптимизирована для обработки реального сигнала так что за раз обрабатывается сразу 2 реальных сигнала которые затем преобразуются в 2 спектра декомпозицией на основе сигналов с четной и нечетной симметрией. Дальнейший синтез единого спектра из двух основан на сдвиге во временной области и соответствующих поворотах в частотной. Комментарии в коде подробно описывают процесс.

spectrum.png

Спасибо за внимание!
The trick.

[Mikle]

Проект при загрузке жалуется на отсутствие comctl32.ocx, хотя такая библиотека в System32 имеется, её регистрация не помогла. Однако, если ошибку игнорировать, то проет загружается и работает, может просто чего-то не хватает? Выглядит как на скриншоте.
На шумовых сигналах весь график спектра хаотично движется вверх-вниз, что выглядит нереально - это бы соответствовало такому же изменению громкости всего спектра частот.

[The trick]

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

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

Проект при загрузке жалуется на отсутствие comctl32.ocx, хотя такая библиотека в System32 имеется, её регистрация не помогла. Однако, если ошибку игнорировать, то проет загружается и работает, может просто чего-то не хватает?

Попробуй поменять версию на такую - Object={6B7E6392-850A-101B-AFC0-4210102A8DA7}#1.0#0; comctl32.Ocx
Она в принципе вообще не нужна, я просто забыл ее убрать оттуда.

[Mikle]

Попробуй поменять версию на такую

Сработало.

Отображение спектра автоматически подстраивается по максимальному уровню сигнала

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

[The trick]

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

Представь белый шум - он имеет равномерную спектральную плотность на всех частотах. Если мы имеем дело с аналоговым сигналом с бесконечным спектром то там амплитуда каждой гармоники стремится к 0, но вклад всех гармоник дает сигнал определенной амплитуды. Для дискретного преобразования в данном случае спектральная плотность равномерно распределяется между отсчетами в частотной области, т.е. мы работаем не с одной гармоникой а группой но все равно - величина спектральной плотности зависит от разрешения. Т.е. если мы берем ДПФ от 16 отсчетов то белый шум равномерно распределится между этими 16 отсчетами и величина каждого будет приблизительно 1/16. Если же мы возьмем ДПФ от 32768 то амплитуда каждого будет 1/32768. Можно подумать что возможно просто нормализовать по количеству отсчетов, но это не подойдет и вот почему.
Сравни 2 сигнала - шум с амплитудой 1 В и синусоиду с амплитудой 1 В. В первом случае гармоники распределены по всей частотной области, во втором только в одном месте. Если делать нормировку по количеству отчетов, то в первом случае мы получим сигнал с очень маленькой амплитудой на каждой частоте, что соответствует истине, но не подходит для визуализации, а синусоида будет соответствовать 0 Дб (если смотреть относительно уровня 1В). Поэтому я выбрал усреднение по пиковому значению, поскольку это показывает именно спектр сигнала независимо от его амплитуды.

[Mikle]

если мы берем ДПФ от 16 отсчетов то белый шум равномерно распределится между этими 16 отсчетами и величина каждого будет приблизительно 1/16

Не совсем. Распределение будет равномерным, если делить по логарифмической шкале, а по Фурье мы получаем деление по линейной шкале (помнится, этот вопрос возникал при написании кругового визуализатора спектра). Или эта проблема решена?
И уровень каждого отсчёта, если их сумма равна 1, будет, по идее, Sqr(1/16), шумы суммируются по такому закону. Речь о действующем, а не амплитудном, значении.

Но в целом мысль понял, действительно, у чистого, долгого синуса уровень будет стремиться к бесконечности, и с этим надо что-то делать. Может, просто давать графику уйти за верхний край окна? Ведь это будет только на очень узких полосах и при чистых, длительных (!) сигналах. Заметь, на синусе с качающейся частотой ширина полосы тем больше, чем быстрее движется опорная частота по шкале. И это не недостаток анализатора, а объективная реальность - у такого сигнала действительно такой спектр.

[The trick]

Не совсем. Распределение будет равномерным, если делить по логарифмической шкале, а по Фурье мы получаем деление по линейной шкале (помнится, этот вопрос возникал при написании кругового визуализатора спектра). Или эта проблема решена?

Так равномерное это и есть линейное, по крайней мере я это и подразумевал. В том анализаторе проблема была с разрешением, в данном анализаторе разрешение по частоте составляет 44100/32768 (максимальное) Гц.

И уровень каждого отсчёта, если их сумма равна 1, будет, по идее, Sqr(1/16), шумы суммируются по такому закону. Речь о действующем, а не амплитудном, значении.

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

Но в целом мысль понял, действительно, у чистого, долгого синуса уровень будет стремиться к бесконечности, и с этим надо что-то делать. Может, просто давать графику уйти за верхний край окна? Ведь это будет только на очень узких полосах и при чистых, длительных (!) сигналах. Заметь, на синусе с качающейся частотой ширина полосы тем больше, чем быстрее движется опорная частота по шкале. И это не недостаток анализатора, а объективная реальность - у такого сигнала действительно такой спектр.

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

[Mikle]

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

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

Прибавь ещё один тип сигнала:

Код: Выделить всё

Private Function Zummer() As Single
  Static d As Single
  Static p As Long

  p = p + 1
  If p > 10000 Then p = p - 20000

  If p > 0 Then
    Zummer = Sin(d)
    d = d + 0.05
    If d > PI Then d = d - PI * 2
  Else
    Zummer = 0
  End If
End Function


У такого сигнала частота не меняется, а спектр меняет ширину.

[The trick]

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

Не пойму про что ты. Смотри. Спектр дискретного сигнала всегда периодический. При увеличении разрешения период будет увеличиваться и так до бесконечности. Когда сигнал бесконечен - это также периодический сигнал только с бесконечным периодом.
Ширина спектральной плотности зависит от того какой сигнал и как он ложится на сетку частот. Когда в реальной задаче мы берем сигнал, то мы берем его на каком-то интервале, а это обозначает что мы умножаем сигнал на функцию прямоугольного окна. Умножение во временной области соответствует свертке в частотной, в итоге вместо единичного пика (дельта функции) мы получаем спектр прямоугольного импульса (sinc), что и показывает твой пример - видно характерный спектр:

sinc.png

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

sinc.gif

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

[Mikle]

Я вообще не про дискретность сигнала. Само собой, что зазубренность графика на твоей картинке - следствие дискретности. Я совершенно о другом, о том, что для определения спектра сигнала нужна его история, представь, что частота дискретизации сигнала неограниченно большая, так же у тебя в фильтре число полос стремится к бесконечности, всё равно при появлении чистого синуса спектр сначала будет широким (рис. 2), потом будет сужаться (рис. 3), а в самый первый момент, который не успеешь заметить, он будет плоским (рис. 1). И все эти спектры не являются суммой спектров синуса и прямоугольного импульса. И это не недостаток метода измерения, это - объективная реальность. Если к передающей антенне подключить в какой-то момент синусоидальный сигнал, то фотоны, излучённые антенной в разные моменты времени, будут иметь частоты в соответствие с этой картиной.

[The trick]

Я вообще не про дискретность сигнала. Само собой, что зазубренность графика на твоей картинке - следствие дискретности. Я совершенно о другом, о том, что для определения спектра сигнала нужна его история, представь, что частота дискретизации сигнала неограниченно большая, так же у тебя в фильтре число полос стремится к бесконечности, всё равно при появлении чистого синуса спектр сначала будет широким (рис. 2), потом будет сужаться (рис. 3), а в самый первый момент, который не успеешь заметить, он будет плоским (рис. 1). И все эти спектры не являются суммой спектров синуса и прямоугольного импульса. И это не недостаток метода измерения, это - объективная реальность. Если к передающей антенне подключить в какой-то момент синусоидальный сигнал, то фотоны, излучённые антенной в разные моменты времени, будут иметь частоты в соответствие с этой эту картиной.

Сигнал о котором ты говоришь уже не будет синусом, т.к. синус определен от минус бесконечности до плюс бесконечности. Во-вторых то о чем ты говоришь уже не является непрерывным спектром, т.к. это оконное преобразование поскольку есть моменты до и после. Т.е. как ты говоришь в первоначальный момент и потом - не применимо к непрерывному преобразованию, т.к. оно определено на отрезке от -∞ до +∞. Т.е. оно включает все моменты времени. Плоский спектр соответствует дельта импульсу - это будет тогда когда только единичный импульс будет в окне. Когда ты говоришь что в первоначальный момент мы будем видеть сначала одно уже обозначает то что ты ограничиваешь сигнал каким-то окном на котором ты измеряешь спектр. А соответственно ограничивание окном есть свертка с функцией окна.

[Mikle]

Сигнал о котором ты говоришь уже не будет синусом, т.к. синус определен от минус бесконечности до плюс бесконечности. Во-вторых то о чем ты говоришь уже не является непрерывным спектром

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

Так равномерное это и есть линейное, по крайней мере я это и подразумевал

Всегда считал, что равномерное - это по логарифмической шкале. Возможно, я ошибаюсь, но тогда получается, что белый шум, который по определению имеет одинаковый уровень спектра на всех частотах, содержит столько мощности в диапазоне от 20 до 40 Гц, сколько и в диапазоне от 15020 до 15040, что выглядит сомнительно.
Ещё по БШ - у тебя там довольно странная формула: (Rnd + Rnd + Rnd + Rnd + Rnd + Rnd) / 3 - 1. Почему не просто (Rnd - 0.5) * 2? Ведь распределение становится неравномерным.

[The trick]

Всегда считал, что равномерное - это по логарифмической шкале. Возможно, я ошибаюсь, но тогда получается, что белый шум, который по определению имеет одинаковый уровень спектра на всех частотах, содержит столько мощности в диапазоне от 20 до 40 Гц, сколько и в диапазоне от 15020 до 15040, что выглядит сомнительно.

Так и есть. Мощность будет одинакова для белого шума на бесконечном промежутке времени.

Ещё по БШ - у тебя там довольно странная формула: (Rnd + Rnd + Rnd + Rnd + Rnd + Rnd) / 3 - 1. Почему не просто (Rnd - 0.5) * 2? Ведь распределение становится неравномерным.

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

[Mikle]

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

Код: Выделить всё

        fValue = m_fSpecBuffer(lIndex)
        '// Get maximum between frequencies
        Do While lPrevIndex <= lIndex

            If m_fSpecBuffer(lPrevIndex) > fValue Then
                fValue = m_fSpecBuffer(lPrevIndex)
            End If

            lPrevIndex = lPrevIndex + 1

        Loop


на эти:

Код: Выделить всё

        If lIndex > lPrevIndex Then
          fValue = m_fSpecBuffer(lIndex)
          For i = lPrevIndex + 1 To lIndex - 1
            fValue = fValue + m_fSpecBuffer(i)
          Next i
          fValue = fValue / (lIndex - lPrevIndex)
          lPrevIndex = lIndex
        Else
          fValue = m_fSpecBuffer(lIndex)
        End If


Стало ровно, хочу ещё убрать этот автоуровень, из-за которого весь график плавает вверх-вниз, пока не разобрался где.
Нашёл баг - посмотри на поведение сигнала Frequency Modulation в разрешении 32768. Не пойму, откуда такое берётся, боюсь проблема в реализации самого FFT, там мне сложнее разобраться.

[The trick]

Стало ровно, хочу ещё убрать этот автоуровень, из-за которого весь график плавает вверх-вниз, пока не разобрался где.

Переменная m_fPeak хранит текущий усредненный пик.

Я заменил в процедуре Refresh эти строки:

Спасибо, посмотрю.

Нашёл баг - посмотри на поведение сигнала Frequency Modulation в разрешении 32768. Не пойму, откуда такое берётся, боюсь проблема в реализации самого FFT, там мне сложнее разобраться.

Не пойму в чем баг. У меня вроде все нормально отображается.

[Mikle]

Не пойму в чем баг. У меня вроде все нормально отображается.

У меня на 32768 и не так явно на 16384 так, на более низких частотах ничего подобного:

[Mikle]

Переменная m_fPeak хранит текущий усредненный пик.

Отлично. Я заменил эту строку:

Код: Выделить всё

fValue = fValue / m_fPeak

на такую:

Код: Выделить всё

fValue = fValue * Sqr(m_lResolution)

Теперь уровень никуда не плывёт, осталось шкалу слева растянуть децибел на 90, и будет то, что надо.

[The trick]

У меня на 32768 и не так явно на 16384 так, на более низких частотах ничего подобного:

Извиняюсь за долгое отсутствие. На самом деле спектр такой и будет, это спектр chirp-сигнала, к тому же нелинейного. К примеру если увеличить частоту модуляции то это станет проявляется и на более низких разрешениях:

Вод даже для сравнения, если повысить частоту до 6 Гц и сравнить мой спектроанализатор со спектроанализатором программы GoldWave:

chirp.png

[Mikle]

Да, я проверил своим спектроанализатором на резонаторах - такой же эффект наблюдается.

Вообще, у меня замысел - сделать спектроанализатор для отстройки музыкальной аппаратуры. Существуют такие приборы в железе, я видел программы под Mac и iOs, хочу сделать под Винду.
Такие системы, сколько я их видел, используют не белый, а розовый шум, при идеальной настройке на экране прямая линия, то есть зависимость чувствительности от частоты не такая, как у твоего анализатора. Можно в твой анализатор ввести поправочные коэффициенты, но спектроанализатор на резонаторах сразу имеет именно такую характеристику, с другой стороны, нагружает процессор сильнее. Пока решаю, на чём остановиться.

[The trick]

У розового шума ведь ослабление на 3db/oct? Прямая линия получится если усреднять на бесконечном интервале. Если можешь, нарисуй что ты имеешь в виду.

[Mikle]

У розового шума ведь ослабление на 3db/oct?

Это по линейной частотной шкале. Выходит, что по логарифмической розовый шум становится равномерным:

Розовый шум характеризуется спектральной плотностью мощности, которая обратно пропорциональна частоте

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

[The trick]

Это по линейной частотной шкале. Выходит, что по логарифмической розовый шум становится равномерным

Ну да, так и есть.

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

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

[Mikle]

Не понятно почему розовый шум.

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

Код: Выделить всё

    Static p As Long

    p = (p + 1) Mod 44100
    fSample = (Sin(200 * p * 2 * PI / 44100) + Sin(600 * p * 2 * PI / 44100) + Sin(2000 * p * 2 * PI / 44100)) / 30


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

[The trick]

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

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

[Mikle]

Я подставлял разные значения частот, картина та же, дело не в попадании между центрами замеряемых полос.
Разобрался, столбы разные потому, что я убрал выравнивание по максимумам, но с этим выравниванием шум выгляди нереалистично, как я выше писал. Ни одно, так другое.
Я так понял, FFT сразу даёт значения по максимумам, а не средневзвешенное? Если так, проблему не решить.
На резонаторах, вроде, с этим нет проблем.
Ещё я не понял, почему размерность массивов вдвое меньше разрешения?

[The trick]

Я подставлял разные значения частот, картина та же, дело не в попадании между центрами замеряемых полос.

Дело в этом. Если частота будет кратна размеру окна, то будут одинаковые пики, т.к. они будут единичные (не распределенные).

Разобрался, столбы разные потому, что я убрал выравнивание по максимумам, но с этим выравниванием шум выгляди нереалистично, как я выше писал. Ни одно, так другое.

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

Я так понял, FFT сразу даёт значения по максимумам, а не средневзвешенное? Если так, проблему не решить.

FFT дает значения спектральной плотности. Что значит по максимумам?

На резонаторах, вроде, с этим нет проблем.

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

Ещё я не понял, почему размерность массивов вдвое меньше разрешения?

Потому что для 32768 используется FFT размерности 16384 которое за один раз вычисляет спектр двух 16384 последовательностей. Полученные спектры разделяются на основе симметрии. Затем полный спектр синтезируется также как это делается в обычном FFT на последнем шаге (дублирование и поворот) из двух "суб-спектров".

[Mikle]

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

Я не разбирался с математикой FFT, поэтому вопрос может быть наивным.
А на FFT этого нет? Если мы подаём чистый синус с частотой конкретного бина, то на соседях это никак не отражается?

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

Я уже как-то показывал эти резонаторы (на круговом фильтре), это численная эмуляция процесса в колебательном контуре (либо груза на пружине, что то же самое). Объём вычислений такого фильтра никак не зависит от добротности, порядок всегда = 1, если под порядком понимать количество последовательно включенных резонаторов.

[The trick]

Если мы подаём чистый синус с частотой конкретного бина, то на соседях это никак не отражается?

Если частота точно равна частоте бина, то да, не отражается (не считая погрешности округления). Если частота сдвинута, то тут уже в зависимости от окна, я кидал выше анимацию.

Я уже как-то показывал эти резонаторы (на круговом фильтре), это численная эмуляция процесса в колебательном контуре (либо груза на пружине, что то же самое).

Нужно посмотреть.