Во многих случаях задача получения (вычисления) спектра сигнала выглядит следующим образом. Имеется АЦП, который с частотой дискретизации Fd преобразует непрерывный сигнал, поступающий на его вход в течение времени Т, в цифровые отсчеты - N штук. Далее массив отсчетов подается в некую программку, которая выдает N/2 каких-то числовых значений (программист, который утянул из инета написал программку, уверяет, что она делает преобразование Фурье).

Чтобы проверить, правильно ли работает программа, сформируем массив отсчетов как сумму двух синусоид sin(10*2*pi*x)+0,5*sin(5*2*pi*x) и подсунем программке. Программа нарисовала следующее:

рис.1 График временной функции сигнала

рис.2 График спектра сигнала

На графике спектра имеется две палки (гармоники) 5 Гц с амплитудой 0.5 В и 10 Гц - с амплитудой 1 В, все как в формуле исходного сигнала. Все отлично, программист молодец! Программа работает правильно.

Это значит, что если мы подадим на вход АЦП реальный сигнал из смеси двух синусоид, то мы получим аналогичный спектр, состоящий из двух гармоник.

Итого, наш реальный измеренный сигнал, длительностью 5 сек , оцифрованный АЦП, то есть представленный дискретными отсчетами, имеет дискретный непериодический спектр.

С математической точки зрения - сколько ошибок в этой фразе?

Теперь начальство решило мы решили, что 5 секунд - это слишком долго, давай измерять сигнал за 0.5 сек.



рис.3 График функции sin(10*2*pi*x)+0,5*sin(5*2*pi*x) на периоде измерения 0.5 сек

рис.4 Спектр функции

Что-то как бы не то! Гармоника 10 Гц рисуется нормально, а вместо палки на 5 Гц появилось несколько каких-то непонятных гармоник. Смотрим в интернетах, что да как…

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

рис.5 Добили нулей до 5 сек

рис.6 Получили спектр

Все равно не то, что было на 5 секундах. Придется разбираться с теорией. Идем в Википедию - источник знаний.

2. Непрерывная функция и представление её рядом Фурье

Математически наш сигнал длительностью T секунд является некоторой функцией f(x), заданной на отрезке {0, T} (X в данном случае - время). Такую функцию всегда можно представить в виде суммы гармонических функций (синусоид или косинусоид) вида:

(1), где:

K - номер тригонометрической функции (номер гармонической составляющей, номер гармоники)
T - отрезок, где функция определена (длительность сигнала)
Ak - амплитуда k-ой гармонической составляющей,
θk- начальная фаза k-ой гармонической составляющей

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

(Более строго, среднеквадратичное отклонение ряда от функции f(x) будет стремиться к нулю, но несмотря на среднеквадратичную сходимость, ряд Фурье функции, вообще говоря, не обязан сходиться к ней поточечно. См. https://ru.wikipedia.org/wiki/Ряд_Фурье .)

Этот ряд может быть также записан в виде:

(2),
где , k-я комплексная амплитуда.

Связь между коэффициентами (1) и (3) выражается следующими формулами:

Отметим, что все эти три представления ряда Фурье совершенно равнозначны. Иногда при работе с рядами Фурье бывает удобнее использовать вместо синусов и косинусов экспоненты мнимого аргумента, то есть использовать преобразование Фурье в комплексной форме. Но нам удобно использовать формулу (1), где ряд Фурье представлен в виде суммы косинусоид с соответствующими амплитудами и фазами. В любом случае неправильно говорить, что результатом преобразования Фурье действительного сигнала будут комплексные амплитуды гармоник. Как правильно говорится в Вики «Преобразование Фурье (ℱ) - операция, сопоставляющая одной функции вещественной переменной другую функцию, также вещественной переменной.»

Итого:
Математической основой спектрального анализа сигналов является преобразование Фурье.

Преобразование Фурье позволяет представить непрерывную функцию f(x) (сигнал), определенную на отрезке {0, T} в виде суммы бесконечного числа (бесконечного ряда) тригонометрических функций (синусоид и\или косинусоид) с определёнными амплитудами и фазами, также рассматриваемых на отрезке {0, T}. Такой ряд называется рядом Фурье.

Отметим еще некоторые моменты, понимание которых требуется для правильного применения преобразования Фурье к анализу сигналов. Если рассмотреть ряд Фурье (сумму синусоид) на всей оси Х, то можно увидеть, что вне отрезка {0, T} функция представленная рядом Фурье будет будет периодически повторять нашу функцию.

Например, на графике рис.7 исходная функция определена на отрезке {-T\2, +T\2}, а ряд Фурье представляет периодическую функцию, определенную на всей оси х.

Это происходит потому, что синусоиды сами являются периодическими функциями, соответственно и их сумма будет периодической функцией.

рис.7 Представление непериодической исходной функции рядом Фурье

Таким образом:

Наша исходная функция - непрерывная, непериодическая, определена на некотором отрезке длиной T.
Спектр этой функции - дискретный, то есть представлен в виде бесконечного ряда гармонических составляющих - ряда Фурье.
По факту, рядом Фурье определяется некоторая периодическая функция, совпадающая с нашей на отрезке {0, T}, но для нас эта периодичность не существенна.

Периоды гармонических составляющих кратны величине отрезка {0, T}, на котором определена исходная функция f(x). Другими словами, периоды гармоник кратны длительности измерения сигнала. Например, период первой гармоники ряда Фурье равен интервалу Т, на котором определена функция f(x). Период второй гармоники ряда Фурье равен интервалу Т/2. И так далее (см. рис. 8).

рис.8 Периоды (частоты) гармонических составляющих ряда Фурье (здесь Т=2π)

Соответственно, частоты гармонических составляющих кратны величине 1/Т. То есть частоты гармонических составляющих Fk равны Fk= к\Т, где к пробегает значения от 0 до ∞, например к=0 F0=0; к=1 F1=1\T; к=2 F2=2\T; к=3 F3=3\T;… Fk= к\Т (при нулевой частоте - постоянная составляющая).

Пусть наша исходная функция, представляет собой сигнал, записанный в течение Т=1 сек. Тогда период первой гармоники будет равен длительности нашего сигнала Т1=Т=1 сек и частота гармоники равна 1 Гц. Период второй гармоники будет равен длительности сигнала, деленной на 2 (Т2=Т/2=0,5 сек) и частота равна 2 Гц. Для третьей гармоники Т3=Т/3 сек и частота равна 3 Гц. И так далее.

Шаг между гармониками в этом случае равен 1 Гц.

Таким образом сигнал длительностью 1 сек можно разложить на гармонические составляющие (получить спектр) с разрешением по частоте 1 Гц.
Чтобы увеличить разрешение в 2 раза до 0,5 Гц - надо увеличить длительность измерения в 2 раза - до 2 сек. Сигнал длительностью 10 сек можно разложить на гармонические составляющие (получить спектр) с разрешением по частоте 0,1 Гц. Других способов увеличить разрешение по частоте нет.

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

3. Дискретные сигналы и дискретное преобразование Фурье

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

Обычная схема измерения и оцифровки сигнала выглядит следующим образом.

рис.9 Схема измерительного канала

Сигнал с измерительного преобразователя поступает на АЦП в течение периода времени Т. Полученные за время Т отсчеты сигнала (выборка) передаются в компьютер и сохраняются в памяти.

рис.10 Оцифрованный сигнал - N отсчетов полученных за время Т

Какие требования выдвигаются к параметрам оцифровки сигнала? Устройство, преобразующее входной аналоговый сигнал в дискретный код (цифровой сигнал) называется аналого-цифровой преобразователь (АЦП, англ. Analog-to-digital converter, ADC) (Wiki).

Одним из основных параметров АЦП является максимальная частота дискретизации (или частота семплирования, англ. sample rate) - частота взятия отсчетов непрерывного во времени сигнала при его дискретизации. Измеряется в герцах. ((Wiki))

Согласно теореме Котельникова, если непрерывный сигнал имеет спектр, ограниченный частотой Fмакс, то он может быть полностью и однозначно восстановлен по его дискретным отсчетам, взятым через интервалы времени , т.е. с частотой Fd ≥ 2*Fмакс, где Fd - частота дискретизации; Fмакс - максимальная частота спектра сигнала. Другими слова частота оцифровки сигнала (частота дискретизации АЦП) должна как минимум в 2 раза превышать максимальную частоту сигнала, который мы хотим измерить.

А что будет, если мы будем брать отсчеты с меньшей частотой, чем требуется по теореме Котельникова?

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

Рис. 11. Появление ложного сигнала низкой частоты при недостаточно высокой частоте дискретизации

Чтобы избежать эффекта алиасинга перед АЦП ставят специальный антиалиасинговый фильтр - ФНЧ (фильтр нижних частот), который пропускает частоты ниже половины частоты дискретизации АЦП, а более высокие частоты зарезает.

Для того, чтобы вычислить спектр сигнала по его дискретным отсчетам используется дискретное преобразование Фурье (ДПФ). Отметим еще раз, что спектр дискретного сигнала «по определению» ограничен частотой Fмакс, меньшей половине частоты дискретизации Fd. Поэтому спектр дискретного сигнала может быть представлен суммой конечного числа гармоник, в отличие от бесконечной суммы для ряда Фурье непрерывного сигнала, спектр которого может быть неограничен. Согласно теореме Котельникова максимальная частота гармоники должна быть такой, чтобы на нее приходилось как минимум два отсчета, поэтому число гармоник равно половине числа отсчетов дискретного сигнала. То есть если в выборке имется N отсчетов, то число гармоник в спектре будет равно N/2.

Рассмотрим теперь дискретное преобразование Фурье (ДПФ).

Сравнивая с рядом Фурье

Видим, что они совпадают, за исключением того, что время в ДПФ имеет дискретный характер и число гармоник ограничено величиной N/2 - половиной числа отсчетов.

Формулы ДПФ записываются в безразмерных целых переменных k, s, где k – номера отсчетов сигнала, s – номера спектральных составляющих.
Величина s показывает количество полных колебаний гармоники на периоде Т (длительности измерения сигнала). Дискретное преобразование Фурье используется для нахождения амплитуд и фаз гармоник численным методом, т.е. «на компьютере»

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

рис.12 Периодическая функция f(x) с периодом Т0, с периодом измерения Т>T0

Как видно на рис.12 функция f(x) периодическая с периодом Т0. Однако из-за того, что длительность измерительной выборки Т не совпадает с периодом функции Т0, функция, получаемая как ряд Фурье, имеет разрыв в точке Т. В результате спектр данной функции будет содержать большое количество высокочастотных гармоник. Если бы длительность измерительной выборки Т совпадала с периодом функции Т0, то в полученном после преобразования Фурье спектре присутствовала бы только первая гармоника (синусоида с периодом равным длительности выборки), поскольку функция f(x) представляет собой синусоиду.

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

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

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

Рис.13 Пример функции и спектра сигнала кинематической погрешности редуктора

При меньшей длительности картина будет выглядеть «хуже»:

Рис.14 Пример функции и спектра сигнала вибрации ротора

На практике бывает сложно понять, где «реальные составляющие», а где «артефакты», вызванные некратностью периодов составляющих и длительности выборки сигнала или «скачками и разрывами» формы сигнала. Конечно слова «реальные составляющие» и «артефакты» не зря взяты в кавычки. Наличие на графике спектра множества гармоник не означает, что наш сигнал в реальности из них «состоит». Это все равно что считать, будто число 7 «состоит» из чисел 3 и 4. Число 7 можно представить в виде суммы чисел 3 и 4 - это правильно.

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

Некоторые итоги

1. Реальный измеренный сигнал, длительностью T сек, оцифрованный АЦП, то есть представленный набором дискретных отсчетов (N штук), имеет дискретный непериодический спектр, представленный набором гармоник (N/2 штук).

2. Сигнал представлен набором действительных значений и его спектр представлен набором действительных значений. Частоты гармоник положительны. То, что математикам бывает удобнее представить спектр в комплексной форме с использованием отрицательных частот не значит, что «так правильно» и «так всегда надо делать».

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

Использованные материалы и другие полезные материалы.

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

Быстрое преобразование Фурье

Быстрое преобразование Фурье - это алгоритм вычисления, который успешно использует свойства периодичности тригонометрических функций для того, чтобы избежать ненужных вычислений в дискретном преобразовании Фурье (ДПФ), за счёт чего быстрее осуществляется поиск коэффициентов в разложении Фурье. Основное отличие от дискретного преобразования заключается лишь в методе вычисления числовых значений (алгоритме), а не в самой обработке сигнала. И в случае БПФ, и в случае ДПФ результат вычислений один и тот же . Единственным требованием для алгоритма БПФ является размер выборки, кратный N = 2L, где L - любое положительное целое число. Наиболее распространёнными алгоритмами БПФ по основанию 2 являются: с прореживанием по времени и с прореживанием по частоте.

В данной работе реализован алгоритм БПФ по основанию 2 с прореживанием по времени (алгоритм Кули - Тьюки). Его легко получить, исследуя некоторые закономерности ДПФ. Введём так называемый поворотный коэффициент:

В этом случае в ДПФ коэффициенты Фурье для ряда значений сигнала {f0,f1,…,fN-1} выражаются соотношением:

Рассмотрим ряд сигнала из 4 значений: {f0,f1,f2,f3}. Представим преобразование Фурье в матричном виде (нормировочный коэффициент 1/N внесён в вектор-столбец Сij в правой части выражения):

Расписав поворотные коэффициенты по формуле Эйлера и определив их значения при k = 0, 1, 2,.. 9, можно построить диаграмму (Рис. 2), из которой видна закономерность повторяющихся коэффициентов.

Рисунок 2. Степенной ряд w для N=4

Подставляя численные значения в (4) получим:

То есть значения w, начиная с w4, равны соотвествующему значению от w0 до w3. Далее перепишем матричное уравнение (4) в нестандартном виде (подобные обозначения введены для наглядности дальнейших операций):

Поменяем местами столбцы матрицы, разделив её на две группы: по чётным f0, f2 и нечётным f1, f3 индексам:

Учтём, что wk+1 = wkw1, тогда выражение (6) перепишется в виде:

Используя соотношения:

Получаем искомые коэффициенты разложения в виде вектора-столбца со значениями ячеек:

Графическое изображение алгоритма (Рис. 3) похоже на бабочку с распахнутыми крыльями, поэтому этот метод вычисления называют «бабочкой».

Рисунок 3. Граф «Бабочка» для ряда из 4 членов

Итак, на первом шаге алгоритма идёт разбиение по чётным и нечётным индексам членов ряда значений сигнала. Затем выполняется граф «бабочка», он состоит из двух ступеней, их количество равно степени двойки объёма выборки (N = 4 = 22). На каждом ступени выполняется по две «бабочки» и их общее количество неизменно. Каждой операции «бабочка» соотвествует одна операция умножения. Для сравнения: в ДПФ с выборкой {f0,f1,f2,f3} операцию умножения необходимо было бы выполнить 4Ч4 = 16 раз, а в случае БПФ всего лишь 4 раза .

Этот ряд может быть также записан в виде:

(2),
где , k-я комплексная амплитуда.

Связь между коэффициентами (1) и (3) выражается следующими формулами:

Отметим, что все эти три представления ряда Фурье совершенно равнозначны. Иногда при работе с рядами Фурье бывает удобнее использовать вместо синусов и косинусов экспоненты мнимого аргумента, то есть использовать преобразование Фурье в комплексной форме. Но нам удобно использовать формулу (1), где ряд Фурье представлен в виде суммы косинусоид с соответствующими амплитудами и фазами. В любом случае неправильно говорить, что результатом преобразования Фурье действительного сигнала будут комплексные амплитуды гармоник. Как правильно говорится в Вики «Преобразование Фурье (?) - операция, сопоставляющая одной функции вещественной переменной другую функцию, также вещественной переменной.»

Итого:
Математической основой спектрального анализа сигналов является преобразование Фурье.

Преобразование Фурье позволяет представить непрерывную функцию f(x) (сигнал), определенную на отрезке {0, T} в виде суммы бесконечного числа (бесконечного ряда) тригонометрических функций (синусоид и\или косинусоид) с определёнными амплитудами и фазами, также рассматриваемых на отрезке {0, T}. Такой ряд называется рядом Фурье.

Отметим еще некоторые моменты, понимание которых требуется для правильного применения преобразования Фурье к анализу сигналов. Если рассмотреть ряд Фурье (сумму синусоид) на всей оси Х, то можно увидеть, что вне отрезка {0, T} функция представленная рядом Фурье будет будет периодически повторять нашу функцию.

Например, на графике рис.7 исходная функция определена на отрезке {-T\2, +T\2}, а ряд Фурье представляет периодическую функцию, определенную на всей оси х.

Это происходит потому, что синусоиды сами являются периодическими функциями, соответственно и их сумма будет периодической функцией.

рис.7 Представление непериодической исходной функции рядом Фурье

Таким образом:

Наша исходная функция - непрерывная, непериодическая, определена на некотором отрезке длиной T.
Спектр этой функции - дискретный, то есть представлен в виде бесконечного ряда гармонических составляющих - ряда Фурье.
По факту, рядом Фурье определяется некоторая периодическая функция, совпадающая с нашей на отрезке {0, T}, но для нас эта периодичность не существенна.

Периоды гармонических составляющих кратны величине отрезка {0, T}, на котором определена исходная функция f(x). Другими словами, периоды гармоник кратны длительности измерения сигнала. Например, период первой гармоники ряда Фурье равен интервалу Т, на котором определена функция f(x). Период второй гармоники ряда Фурье равен интервалу Т/2. И так далее (см. рис. 8).

рис.8 Периоды (частоты) гармонических составляющих ряда Фурье (здесь Т=2?)

Соответственно, частоты гармонических составляющих кратны величине 1/Т. То есть частоты гармонических составляющих Fk равны Fk= к\Т, где к пробегает значения от 0 до?, например к=0 F0=0; к=1 F1=1\T; к=2 F2=2\T; к=3 F3=3\T;… Fk= к\Т (при нулевой частоте - постоянная составляющая).

Пусть наша исходная функция, представляет собой сигнал, записанный в течение Т=1 сек. Тогда период первой гармоники будет равен длительности нашего сигнала Т1=Т=1 сек и частота гармоники равна 1 Гц. Период второй гармоники будет равен длительности сигнала, деленной на 2 (Т2=Т/2=0,5 сек) и частота равна 2 Гц. Для третьей гармоники Т3=Т/3 сек и частота равна 3 Гц. И так далее.

Шаг между гармониками в этом случае равен 1 Гц.

Таким образом сигнал длительностью 1 сек можно разложить на гармонические составляющие (получить спектр) с разрешением по частоте 1 Гц.
Чтобы увеличить разрешение в 2 раза до 0,5 Гц - надо увеличить длительность измерения в 2 раза - до 2 сек. Сигнал длительностью 10 сек можно разложить на гармонические составляющие (получить спектр) с разрешением по частоте 0,1 Гц. Других способов увеличить разрешение по частоте нет.

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

3. Дискретные сигналы и дискретное преобразование Фурье

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

Обычная схема измерения и оцифровки сигнала выглядит следующим образом.

рис.9 Схема измерительного канала

Сигнал с измерительного преобразователя поступает на АЦП в течение периода времени Т. Полученные за время Т отсчеты сигнала (выборка) передаются в компьютер и сохраняются в памяти.

рис.10 Оцифрованный сигнал - N отсчетов полученных за время Т

Какие требования выдвигаются к параметрам оцифровки сигнала? Устройство, преобразующее входной аналоговый сигнал в дискретный код (цифровой сигнал) называется аналого-цифровой преобразователь (АЦП, англ. Analog-to-digital converter, ADC) (Wiki).

Одним из основных параметров АЦП является максимальная частота дискретизации (или частота семплирования, англ. sample rate) - частота взятия отсчетов непрерывного во времени сигнала при его дискретизации. Измеряется в герцах. ((Wiki))

Согласно теореме Котельникова, если непрерывный сигнал имеет спектр, ограниченный частотой Fмакс, то он может быть полностью и однозначно восстановлен по его дискретным отсчетам, взятым через интервалы времени , т.е. с частотой Fd ? 2*Fмакс, где Fd - частота дискретизации; Fмакс - максимальная частота спектра сигнала. Другими слова частота оцифровки сигнала (частота дискретизации АЦП) должна как минимум в 2 раза превышать максимальную частоту сигнала, который мы хотим измерить.

А что будет, если мы будем брать отсчеты с меньшей частотой, чем требуется по теореме Котельникова?

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

Рис. 11. Появление ложного сигнала низкой частоты при недостаточно высокой частоте дискретизации

Чтобы избежать эффекта алиасинга перед АЦП ставят специальный антиалиасинговый фильтр - ФНЧ (фильтр нижних частот), который пропускает частоты ниже половины частоты дискретизации АЦП, а более высокие частоты зарезает.

Для того, чтобы вычислить спектр сигнала по его дискретным отсчетам используется дискретное преобразование Фурье (ДПФ). Отметим еще раз, что спектр дискретного сигнала «по определению» ограничен частотой Fмакс, меньшей половине частоты дискретизации Fd. Поэтому спектр дискретного сигнала может быть представлен суммой конечного числа гармоник, в отличие от бесконечной суммы для ряда Фурье непрерывного сигнала, спектр которого может быть неограничен. Согласно теореме Котельникова максимальная частота гармоники должна быть такой, чтобы на нее приходилось как минимум два отсчета, поэтому число гармоник равно половине числа отсчетов дискретного сигнала. То есть если в выборке имется N отсчетов, то число гармоник в спектре будет равно N/2.

Рассмотрим теперь дискретное преобразование Фурье (ДПФ).

Сравнивая с рядом Фурье

Видим, что они совпадают, за исключением того, что время в ДПФ имеет дискретный характер и число гармоник ограничено величиной N/2 - половиной числа отсчетов.

Формулы ДПФ записываются в безразмерных целых переменных k, s, где k – номера отсчетов сигнала, s – номера спектральных составляющих.
Величина s показывает количество полных колебаний гармоники на периоде Т (длительности измерения сигнала). Дискретное преобразование Фурье используется для нахождения амплитуд и фаз гармоник численным методом, т.е. «на компьютере»

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

рис.12 Периодическая функция f(x) с периодом Т0, с периодом измерения Т>T0

Как видно на рис.12 функция f(x) периодическая с периодом Т0. Однако из-за того, что длительность измерительной выборки Т не совпадает с периодом функции Т0, функция, получаемая как ряд Фурье, имеет разрыв в точке Т. В результате спектр данной функции будет содержать большое количество высокочастотных гармоник. Если бы длительность измерительной выборки Т совпадала с периодом функции Т0, то в полученном после преобразования Фурье спектре присутствовала бы только первая гармоника (синусоида с периодом равным длительности выборки), поскольку функция f(x) представляет собой синусоиду.

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

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

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

Рис.13 Пример функции и спектра сигнала кинематической погрешности редуктора

При меньшей длительности картина будет выглядеть «хуже»:

Рис.14 Пример функции и спектра сигнала вибрации ротора

На практике бывает сложно понять, где «реальные составляющие», а где «артефакты», вызванные некратностью периодов составляющих и длительности выборки сигнала или «скачками и разрывами» формы сигнала. Конечно слова «реальные составляющие» и «артефакты» не зря взяты в кавычки. Наличие на графике спектра множества гармоник не означает, что наш сигнал в реальности из них «состоит». Это все равно что считать, будто число 7 «состоит» из чисел 3 и 4. Число 7 можно представить в виде суммы чисел 3 и 4 - это правильно.

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

Некоторые итоги

1. Реальный измеренный сигнал, длительностью T сек, оцифрованный АЦП, то есть представленный набором дискретных отсчетов (N штук), имеет дискретный непериодический спектр, представленный набором гармоник (N/2 штук).

2. Сигнал представлен набором действительных значений и его спектр представлен набором действительных значений. Частоты гармоник положительны. То, что математикам бывает удобнее представить спектр в комплексной форме с использованием отрицательных частот не значит, что «так правильно» и «так всегда надо делать».

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

Использованные материалы и другие полезные материалы.

В статье описан небольшой анализатор аудиоспектра (0 - 10 кГц), состоящий из ЖК-дисплея 16x2 и микроконтроллера ATmega32. Используется простой алгоритм ДПФ (Дискретное Преобразование Фурье). БПФ (Быстрое Преобразование Фурье) отличается от ДПФ только большей скоростью но и более сложным алгоритмом.

ДПФ медленный по сравнению с БПФ. Данный ЖК анализатор спектра не требует большой скорости, которую может обеспечить БПФ, и если изображение на экране будет меняться с частотой около 30 кадров / сек, то это более чем достаточно для визуализации звукового спектра. Для ЖК-дисплея не рекомендуется слишком высокая частота обновления. Звук с частотой дискретизации 20 кГц даёт 32 точки ДПФ. Поскольку результат преобразования симметричен, мне нужно использовать только первые 16 результатов. Соответственно максимальная частота 10 кГц. Таким образом, 10кГц/16 = 625Гц.

Можно увеличить скорость вычисления ДПФ. Если есть точка N ДПФ, то необходимо найти синус и косинус (N ^ 2) / 2. Для 32-точечного ДПФ, необходимо найти синус и косинус 512. Прежде чем искать синус и косинус, нам нужно найти угол (градусы), который занимает некоторое время процессора. Для ускорения этого сделаны таблицы для синуса и косинуса. Синус и косинус сделаны 16-битными переменными, умножив значения синуса и косинуса на 10000. После преобразования нужно разделить каждый результат на 10000. Теперь можно рассчитать 120 32-точечных ДПФ в секунду, что более чем достаточно для анализатора спектра.

Дисплей

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

Аудио вход

Одной из наиболее важных частей анализатора спектра является получение сигнала с электретного микрофона. Особое внимание должно быть уделено разработке предварительного усилителя для микрофона. Нам нужно установить нулевой уровень на входе АЦП и максимальный уровень равный половине напряжения питания, т.е. 2,5В. На него может подаваться напряжение от -2,5В до +2,5В. Предусилитель должен быть настроен так, чтобы не превышать этих границ. В схеме использован операционный усилитель LM324 в качестве предварительного усилителя для микрофона.

В качестве устройства отображения используется двухстрочный символьный ЖК индикатор. Основным моментом при реализации данного проекта является не аппаратная часть, а программная, точнее реализация дискретного преобразования Фурье (ДПФ) на 8-разрядном микроконтроллере. Сразу следует отметить, что автор не является экспертом в этой области и поэтому начал с основ - с простого дискретного преобразования Фурье. Алгоритм быстрого преобразования Фурье является не только быстрым, но и достаточно сложным.

Дискретное преобразование Фурье (в англоязычной литературе DFT, Discrete Fourier Transform) - это одно из преобразований Фурье, широко применяемых в алгоритмах цифровой обработки сигналов (его модификации применяются в сжатии звука в MP3, сжатии изображений в JPEG и др.), а также в других областях, связанных с анализом частот в дискретном (к примеру, оцифрованном аналоговом) сигнале. Дискретное преобразование Фурье требует в качестве входа дискретную функцию. Такие функции часто создаются путем дискретизации (выборки значений из непрерывных функций).

Принципиальная схема анализатора спектра звукового сигнала очень проста и условно ее можно разделить на цифровую часть и аналоговую.

Цифровая часть образована микроконтроллером и подключенным к нему ЖК индикатором. Микроконтроллер тактируется от кварцевого резонатора 16 МГц, в качестве опорного напряжения для АЦП микроконтроллера используется напряжение питания +5 В.
Шина данных ЖК индикатора подключена к порту C микроконтроллера (линии ввода/вывода PC0-PC3), шина управления подключена к порту D(PD5, PD6) микроконтроллера. Индикатор работает в 4-битном режиме. Переменный резистор номиналом 4.7 кОм используется для регулировки контрастности. Для работы с индикатором были созданы пользовательские символы для отображения 8 горизонтальных столбиков анализатора, эти пользовательские символы занимают все 64 Байта ОЗУ ЖК индикатора.

Микроконтроллер работает от внешнего кварцевого резонатора 16 МГц.

Аналоговая часть устройства - это самая важная часть и представляет собой предварительный усилитель сигнала электретного микрофона, выход которого подключается к каналу ADC0 встроенного в микроконтроллер АЦП. Уровень нуля на входе АЦП нам необходимо установить равным точно половине опорного напряжения, т.е. 2.5 В. В этом случае мы сможем использовать положительную и отрицательную полуволну сигнала, но его амплитуда не должна превышать установленный предел, т.е. коэффициент усиления должен быть точно настроен для предотвращения перегрузки. Всем вышеуказанным условиям отвечает распространенная микросхема низкопотребляющего операционного усилителя .

Алгоритм ДПФ несколько медленнее в сравнении с быстрым преобразованием Фурье. Но наш анализатор спектра не требует высокой скорости, и если он способен обеспечить скорость обновления около 30 кадров в секунду, этого будет более чем достаточно для визуализации спектра звукового сигнала. В любом случае, в нашем варианте возможно достичь скорости 100 кадров в секунду, но это уже слишком высокое значение параметра для двухстрочного символьного ЖК индикатора и оно не рекомендуется. Частота дискретизации равна 20 кГц для 32 точечного дискретного преобразования Фурье и поскольку результат преобразования симметричен, нам нужно использовать только первую половину, т.е. первые 16 результатов. Следовательно, мы можем отображать частотный спектр в диапазоне до 10 кГц и разрешение анализатора составляет 10 кГц/16 = 625 Гц.

Автором конструкции были предприняты попытки увеличения скорости вычисления ДПФ. Если это преобразование имеет N точек, то мы должны найти N2/2 значений синуса и косинуса. Для нашего 32 точечного преобразования необходимо найти 512 значений синуса и косинуса. Но, прежде чем найти их нам необходимо вычислить угол (градусы), что займет некторое процессорное время, поэтому было решено использовать для этих вычислений таблицы значений. При расчетах в программе микроконтроллера не используются вычисления с плавающей точкой и числа двойной точности (double), так как это займет больше времени на обработку на 8-разрядном микроконтроллере. Вместо этого значения в таблицах поиска используются 16-разрядные данные целочисленного типа (integer), умноженные на 10000. Затем после выполнения преобразования результаты делятся на 10000. При таком подходе имеется возможность выполнять 120 32-точечных преобразований в секунду, что более чем достаточно для нашего устройства.

Демонстрация работы анализатора спектра на микроконтроллере ATmega32

Загрузки

Исходный код (программа микроконтроллера, таблицы данных синуса, косинуса и угла) -

• Понятно, что на АВР-ке дальше светомузыки сложно уехать, не те параметры. Но 120 32-точечных преобразований в секунду для большинства задач может быть достаточно. А выборку 625Гц, можно конечно и другую взять, по точнее потеряв частоту обновления. Стоит отметить, что МК будет себя плохо чувствовать, в плане производительности мало что на него еще навешаешь. Но тут можно же организовать выдачу результата по аппаратным методам передачи данных. Тогда это будет вспомогательный МК, а основной будет только принимать с него данный и обрабатывать совместимо с другими процессами. По большому счету все же в частоту проца упирается. Когда-то получалось разгонять мегу выше 20 МГц, но для этих задач наверно получим только глюки на высоких частотах. Идея хороша, только бы больше мат части расписано было бы... именно ее реализация на МК
• я и поинтересней анализаторы делал: You Tube или вариант на цветном ЖКИ: You Tube в основе знаменитая библиотека Чена:)
• "нам необходимо вычислить угол (градусы)" А может кто-нибудь подробнее рассказать как рассчитываются значения для этих таблиц?
• С таблицей синусов и косинусов все понятно. Не понятно как рассчитываются значения в таблице degree_lookup?