По работе неоднократно сталкивался с необходимостью быстро определить наличие в сигнале гармонических составляющих. Часто для примерной оценки достаточно воспользоваться алгоритмом быстрого преобразования Фурье. Тем более, что его реализации есть практически во всех математических пакетах и библиотеках, да и собственноручно реализовать не составит особого труда. Между тем, опыт показывает, что, при всей своей простоте, метод начинает вызывать некоторые вопросы, когда возникает необходимость не просто посмотреть наличие дискреток в сигнале, но и выяснить их абсолютные значения, т.е. нормализовать полученный результат.
В этой статье я постараюсь объяснить, что же все-таки выдает в качестве результата fft (Fast Fourier transform) на примере MATLAB (и в качестве бонуса проведу небольшой ликбез по этому весьма полезному, на мой взгляд, языку).
MATLAB позволяет не заморачиваться с ручным удалением ненужных объектов, однако, при работе с более менее объемными массивами данных, имеет привычку капризничать и жаловаться на недостаток памяти. Для освобождения памяти используется процедура clear с указанием имени объекта, который необходимо удалить.
AGalilov: Преобразование Фурье «на пальцах»
С этого и начнем. Так как все необходимое мы сгенерируем самостоятельно, можно смело удалять все, что накопилось в рабочем пространстве за активную сессию, просто добавив ключевое слово all:
clear all% Очистка памяти
Итак, прежде всего, зададим исходные данные для нашей модели. Фурье анализ идеально подходит для выделения гармонических сигналов на фоне помех. Для того чтобы продемонстрировать это, возьмем в качестве сигнала сумму некоторой постоянной и двух синусоид с разной частотой и амплитудой. Дисперсию шума возьмем в 3 раза больше амплитуды первой синусоиды.
Так же зададим количество частотных полос, которые должен будет посчитать fft алгоритм. Точка с запятой в конце строк не обязательна, и если ее не ставить, результат вычисления функций и задания переменных будет дублироваться в командную строку, что можно использовать для отладки кода (однако, 512 значений сплошным полотном в командной строке вряд ли помогут вам, тем более что их вывод тоже занимает некоторое количество времени, так что все же лучше не забывать закрывать строки).
%% Параметры
Tm=5;% Длина сигнала (с)
Fd=512;% Частота дискретизации (Гц)
Ak=0.5;% Постоянная составляющая (Попугаев)
A1=1;% Амплитуда первой синусоиды (Попугаев)
A2=0.7;% Амплитуда второй синусоиды (Попугаев)
F1=13;% Частота первой синусоиды (Гц)
F2=42;% Частота второй синусоиды (Гц)
Phi1=0;% Начальная фаза первой синусоиды (Градусов)
Phi2=37;% Начальная фаза второй синусоиды (Градусов)
An=3*A1;% Дисперсия шума (Попугаев)
FftL=1024;% Количество линий Фурье спектра
MATLAB (Matrix Laboratory), как следует из названия, предназначен прежде всего для работы с массивами, практически все алгоритмы счета в нем оптимизированы для работы с векторами. Обилие удобных инструментов работы так же ненавязчиво подталкивает представлять как можно больше исходных данных в виде матриц. В частности, можно легко сгенерировать массив возрастающих (убывающих) величин с заданным шагом (1Fd в данном примере):
Основы ЦОС: 18. Преобразование Фурье (ссылки на скачивание скриптов в описании)
%% Генерация рабочих массивов
T=0:1/Fd:Tm;% Массив отсчетов времени
Случайный Гауссов шум задается функцией randn, результатом которой является массив размерности, заданной в ее параметрах. Для единообразия зададим его в виде строки (первый параметр 1) длиной соответствующей длине нашего массива отсчетов времен, что в свою очередь вычислим функцией length.
Noise=An*randn(1,length(T));% Массив случайного шума длиной равной массиву времени
Символ * используется для обозначения перемножения. Т.к. чаще всего действия производятся над векторами, то и умножение подразумевается векторное, но так же легко можно, не перегружая этот оператор, использовать его для поэлементного перемножения, добавив перед ним точку (.*). При умножении вектора на скаляр точка перед умножением не является обязательной. Также добавленная точка сделает поэлементным возведение в степень и деление.
Signal=Ak+A1*sind((F1*360).*T+Phi1)+A2*sind((F2*360).*T+Phi2);% Массив сигнала (смесь 2х синусоид и постоянной составляющей)
Теперь перейдем к тому, ради чего и затевалась данная статья — функции fft(). Аргументами стандартной функции MATLAB являются сам сигнал (в нашем случае Signal), размерность вектора-результата (FftL), а также измерение.
Последний аргумент определяет вдоль какого измерения располагается сигнал в случае если на вход подается многомерный массив (Иногда этот параметр ошибочно принимают за размерность преобразования Фурье, но это не так. Хотя в MATLAB есть реализации 2-хмерного fft2() и многомерного fftn() алгоритмов). Так как наш сигнал представляет собой вектор, то его вполне можно опустить.
Рассмотрим сначала сигнал без примеси шума. В качестве результата мы получим вектор комплексных чисел. Это и есть представление нашего сигнала в частотном домене в показательной форме. Т.е. модули этих комплексных чисел представляют амплитуды соответствующих частот (точнее полосы частот см. дальше), а аргументы – их начальные фазы. И если полученная фаза, однозначно вычисляется в радианах, то с амплитудой и частотами не все так просто.
Например, если мы просто применим к нашему сигналу преобразование Фурье, и возьмем абсолютные значения вектора на выходе, то получим приблизительно следующую картинку:
Для построения двухмерных графиков удобно использовать функцию plot. Основные параметры, используемые в этой функции – одномерные массивы точек, первый задает ось ординат, второй – значение функции в соответствующих точках. Если передать только один массив, то он будет отображен с фиксированным шагом 1.
Если присмотреться к полученной картинке выяснится, что она несколько отличается от наших ожиданий. На приведенном графике 5 пиков вместо ожидаемых 3х (постоянная + 2 синусоиды), их амплитуды не совпадают с амплитудами исходных сигналов, и ось абсцисс вряд ли отображает частоты.
Прежде всего, следует учитывать, что счет алгоритма устроен таким образом, что перебираются не только положительные, но и отрицательные частоты и правая часть графика является «зеркальным» отображением реального спектра. Т.е. на самом деле 0 (которому соответствует постоянная часть сигнала) должен приходиться на середину массива. Ситуацию можно поправить, совершив циклический сдвиг на половину длины массива. Для этих целей, в MATLAB существует функция сдвига fftshift(), смещающая первый элемент в середину массива:
Теперь обратим внимание на ось значений.
Согласно теореме отсчетов (так же известной как теорема Найквиста-Шеннона или более патриотично теорема Котельникова) спектр дискретного сигнала будет ограничен половиной частоты дискретизации (Fd). Или в нашем случае –Fd/2 слева и Fd/2 справа. Т.е. весь полученный массив покрывает Fd частот. Отсюда, учитывая что мы знаем (вернее даже самостоятельно задали в качестве параметра) длину массива, получим частоты в виде массива значений от –Fd/2 до Fd/2 с шагом Fd/FftL (на самом деле крайняя правая частота будет меньше границы на один отсчет т.е. Fd/2-Fd/FftL):
Если посмотреть на фазы частот, можно заметить, что они равны отрицательным фазам соответствующих отрицательных частот. Учитывая равенство амплитуд левой и правой частей спектра и соответствие их фаз с точностью до знака, весь спектр будет эквивалентен своей положительной части с удвоенной амплитудой. Исключение составляет только 0 элемент, который не имеет зеркальной половины. Таким образом, можно избавиться от «непонятных» и зачастую ненужных отрицательных частот. Это можно было сделать сразу просто отбросив конец исходного массива и помножив оставшиеся элементы на 2 (за исключением постоянной составляющей):
Теперь это уже похоже на тот результат, который мы ожидаем. Единственное, что смущает теперь – это амплитуды. С этим все достаточно просто. Т.к. быстрое преобразование Фурье фактически представляет собой суммирование сигнала перемноженного на ядро преобразования (комплексную экспоненту) для каждой из частот, то реальный результат будет меньше полученного ровно в количество суммирований (частот в результате), т.е. полученный результат надо разделить на количество элементов в результате (не забываем, что под результатом понимается весь ответ, вместе с отброшенной частью, т.е. наше заданное FftL):
Стоит упомянуть еще одну вещь. В спектральном представлении вычисляется не значение сигнала на той частоте на которую попал алгоритм (как мы помним частоты следуют с шагом Fd/FftL), а значение в полосе (шириной равной шагу). Т.е. если в эту полосу попало несколько дискреток, то они суммируются. Для примера можно уменьшить количество линий в результате работы алгоритма:
Однако это не означает, что стоит сразу бездумно увеличивать точность работы, т.к. это тоже приводит к негативным последствиям, т.к. если разрешение будет сопоставимо с частотой дискретизации сигнала, в спектр полезут гармоники «окна», которые имеют отношение не к реальному сигналу, а к его дискретному представлению:
Или более близко окрестности одной из дискреток:
Код для нормировки fft будет выглядеть приблизительно следующим образом:
%% Спектральное представление сигнала
FftS=abs(fft(Signal,FftL));% Амплитуды преобразования Фурье сигнала
FftS=2*FftS./FftL;% Нормировка спектра по амплитуде
FftS(1)=FftS(1)/2;% Нормировка постоянной составляющей в спектре
FftSh=abs(fft(Signal+Noise,FftL));% Амплитуды преобразования Фурье смеси сигнал+шум
FftSh=2*FftSh./FftL;% Нормировка спектра по амплитуде
FftSh(1)=FftSh(1)/2;% Нормировка постоянной составляющей в спектре
Нам осталось только вывести результаты. Функция subplot позволяет разбить окно на несколько областей для отображения графиков.
%% Построение графиков
subplot(2,1,1);% Выбор области окна для построения
plot(T,Signal);% Построение сигнала
title(‘Сигнал’);% Подпись графика
xlabel(‘Время (с)’);% Подпись оси х графика
ylabel(‘Амплитуда (Попугаи)’);% Подпись оси у графика
subplot(2,1,2);% Выбор области окна для построения
plot(T,Signal+Noise);% Построение смеси сигнал+шум
title(‘Сигнал+шум’);% Подпись графика
xlabel(‘Время (с)’);% Подпись оси х графика
ylabel(‘Амплитуда (Попугаи)’);% Подпись оси у графика
F=0:Fd/FftL:Fd/2-1/FftL;% Массив частот вычисляемого спектра Фурье
figure% Создаем новое окно
subplot(2,1,1);% Выбор области окна для построения
plot(F,FftS(1:length(F)));% Построение спектра Фурье сигнала
title(‘Спектр сигнала’);% Подпись графика
xlabel(‘Частота (Гц)’);% Подпись оси х графика
ylabel(‘Амплитуда (Попугаи)’);% Подпись оси у графика
subplot(2,1,2);% Выбор области окна для построения
plot(F,FftSh(1:length(F)));% Построение спектра Фурье сигнала
title(‘Спектр сигнала’);% Подпись графика
xlabel(‘Частота (Гц)’);% Подпись оси х графика
ylabel(‘Амплитуда (Попугаи)’);% Подпись оси у графика
Результат выполнения кода будет выглядеть следующим образом:
Несмотря на то, что полезного сигнала не видно на фоне шума, спектральная характеристика позволяет определить его частоту и амплитуду.
Надеюсь данный текст был вам полезен.
- MATLAB
- fft
- быстрое преобразование Фурье
Источник: habr.com
PARALLEL.RU — Информационно-аналитический центр по параллельным вычислениям
Библиотека FFTW для выполнения быстрых преобразований Фурье (БПФ)
- Сайт разработчиков библиотеки — www.fftw.org.
- Документация к версии 2.1.4 библиотеки.
Введение
- Одномерное преобразование Фурье для комплексных чисел
- Многомерное преобразование Фурье для комплексных чисел
- Одномерное преобразование Фурье для действительных чисел
- Многомерное преобразование Фурье для действительных чисел
Работа с FFTW
Сначала происходит построение «плана», который оптимизирует время вычислений для данной задачи. Затем построенный план передается в качестве параметра функциям, которые непосредственно отвечают за вычисление БПФ.
Одномерное преобразование Фурье для комплексных чисел:
#include #define N 100 int main(void) < int i; fftw_complex in[N], out[N]; fftw_plan plan,plan_inv; /* создать оценочный план для прямого БПФ на массиве из N точек */ plan = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE); /* создать оценочный план для обратного БПФ на массиве из N точек */ plan_inv = fftw_create_plan(N, FFTW_BACKWARD, FFTW_ESTIMATE); /* вычислить спектр сигнала, заданного массивом in, в соответствии с планом plan */ fftw_one(plan, in, out); /* вычислительный код программы . */ /* вычислить сигнала по спектру, заданному массивом out, в соответствии с планом plan_inv */ fftw_one(plan_inv, out, in); /* после обратного преобразования Фурье, полученные данные следует нормировать на длину исходного массива */ for(i=0; i < N; i++) < in[i] /= N; >/* удаление ранее созданного плана */ fftw_destroy_plan(plan); fftw_destroy_plan(plan_inv); return 0; >
Пояснения к функции fftw_create_plan
fftw_create_plan(int size, fftw_direction dir, int flags)
- size — количество точек, используемых для вычисления БПФ.
- dir — флаг прямого FFTW_FORWARD или обратного FFTW_BACKWARD преобразования Фурье.
- flag — FFTW_ESTIMATE создает план, который не является полностью оптимальным для данной конфигурации; FFTW_MEASURE — используется для создания оптимального плана, в этом случае значительно увеличивается время работы по вычислению плана и БПФ.
Также можно использовать флаг FFTW_IN_PLACE:
plan = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_IN_PLACE)
тогда в качестве out указывается значение NULL:
fftw_one(plan, in, NULL)
Компиляция под Unix
Как правило, для сборки исполняемого модуля нужно указать а) директорию, где находятся заголовочные файлы (параметр -I); б) директорию, где находятся библиотечные файлы (параметр -L); в) один или несколько библиотечных файлов (параметр -l); Например:
mpicc -c -I/usr/local/fftw/include offtw.c mpicc offtw.o -L/usr/local/fftw/lib -lfftw -lm
Многомерные преобразования для комплексных чисел
#include #define N 100 #define M 100 #define K 100 int main(void) < int i; fftw_complex *in; fftwnd_plan plan,plan_inv; in=(fftw_complex*)malloc((K*(N*M+M)+K)*sizeof(fftw_complex)); /* создать оценочный план для прямого БПФ на массиве из N точек */ plan = fftw3d_create_plan(N, M, K, FFTW_FORWARD, FFTW_ESTIMATE | FFTW_IN_PLACE); /* создать оценочный план для обратного БПФ на массиве из N точек */ plan_inv = fftw3d_create_plan(N, M, K, FFTW_BACKWARD, FFTW_ESTIMATE | FFTW_IN_PLACE); /* вычислить спектр сигнала, заданного массивом in, в соответствии с планом plan */ fftwnd_one(plan, in, NULL); /* вычислительный код программы . */ /* вычислить сигнала по спектру, заданному массивом in, в соответствии с планом plan_inv */ fftwnd_one(plan_inv, in, NULL); /* после обратного преобразования Фурье, полученные данные следует нормировать на длину исходной области */ for(i=0; i < K*(N*M+M)+K; i++) < in[i].re/= N*M*K; in[i].im/= N*M*K; >/* удаление ранее созданного плана */ fftwnd_destroy_plan(plan); fftwnd_destroy_plan(plan_inv); return 0; >
- fftwnd_plan fftw2d_create_plan(int nx, int ny, fftw_direction dir, int flags) — для двумерного БПФ;
- fftwnd_plan fftw3d_create_plan(int nx, int ny, int nz, fftw_direction dir, int flags) — для трехмерного БПФ.
Здесь nx, ny, nz — количество точек по соответствующим декартовым координатам (x, y, z); Флаги dir и flag имеют те же значения, что и в случае одномерного БПФ.
Функция преобразования сигнал/спектр:
void fftwnd_one(fftwnd_plan plan, fftw_complex *in, fftw_complex *out);
Функция удаления созданного плана:
fftwnd_destroy_plan(fftwnd_plan plan)
Преобразование Фурье для действительных чисел
В FFTW входит набор функций для выполнения БПФ с действительными числами. При сборке программы нужно добавить специальную библиотеку с помощью флага -lrfftw, в дополнение к общей библиотеке FFTW (флаг -lfftw). Например:
mpicc -I/usr/local/fftw/include rfftw_test.c -L/usr/local/fftw/lib -lrfftw -lfftw -lm
Модуль параллельной обработки БПФ для MPI
Этот модуль позволяет работать на многопроцессорных машинах с общей и распределенной памятью. Основное его отличие от однопроцессорного модуля FFTW заключается в том, что на каждом процессоре обрабатывается свое подмножество точек.
Функция создания плана для параллельной обработки принимает, в дополнение к стандартным параметрам, значение коммуникатора для набора процессов (MPI_COMM_WORLD или другой коммуникатор). Обмен данными производиться с помощью функций MPI_Alltoall или MPI_Alltoallv в зависимости от алгоритма распределения данных по процессорам. Надо отметить, что наилучшее распределение по процессорам происходит при условии, что число точек кратно P 2 , где P — число процессов.
В исходный текст программы нужно включить заголовочный файл в дополнение к . При сборке программы нужно добавить специальную библиотеку с помощью флага -lfftw_mpi, в дополнение к общей библиотеке FFTW (флаг -lfftw).
Пример выполнения одномерного комплексного БПФ с помощью параллельной версии FFTW:
#include #include #include #include #include int main(int argc,char **argv) < int size,k,i,rank,Nx=1000; int local_n,local_start,local_n_after_transform,local_start_after_transform,total_local_size; double time; fftw_complex *A; fftw_mpi_plan plan,pinv; MPI_Init(argv); MPI_Comm_size(MPI_COMM_WORLD, MPI_Comm_rank(MPI_COMM_WORLD, /* создать план для прямого БПФ, в MPI версии FFTW указывать флаг FFTW_IN_PLACE не нужно*/ plan = fftw_mpi_create_plan(MPI_COMM_WORLD, Nx, FFTW_FORWARD, FFTW_ESTIMATE); /* функция определяющая длину и порядок данных на конкретном процессором */ fftw_mpi_local_sizes(plan, local_start, local_start_after_transform, A=(fftw_complex*)malloc(total_local_size*sizeof(fftw_complex)); /* обнуление массива :. */ /* заполнение массива данными для дальнейшего счета */ for(k=local_start,i=0;k /* выполнение прямого БПФ */ fftw_mpi(plan,1,A,NULL); pinv = fftw_mpi_create_plan(MPI_COMM_WORLD, Nx, FFTW_BACKWARD, FFTW_ESTIMATE); fftw_mpi_local_sizes(pinv, local_start, local_start_aftеr_transform, fftw_mpi(pinv,1,A,NULL); for(i=0;i < total_local_size;i++) < A[i].re/=Nx; A[i].im/=Nx; >fftw_mpi_destroy_plan(plan). fftw_mpi_destroy_plan(pinv). MPI_Finalize(); return 0; >
Комментарии к используемым функциям
Функция создания плана для одномерного комплексного БПФ:
fftw_mpi_plan fftw_mpi_create_plan(MPI_Comm comm, int nx, fftw_direction dir, int flags);
- comm — коммуникатор группы процессоров, на которых будет производиться вычисление БПФ;
- nx — длина массива данных на всех процессорах;
- dir — флаг прямого FFTW_FORWARD или обратного FFTW_BACKWARD преобразования Фурье.
- flag — FFTW_ESTIMATE создает план, который не является полностью оптимальным для данной конфигурации; FFTW_MEASURE — используется для создания оптимального плана, в этом случае значительно увеличивается время работы по вычислению плана и БПФ. FFTW_IN_PLACE не указывается, т.к. в этом случае он стоит по умолчанию. Если не важен порядок выходных данных можно использовать FFTW_SCRAMBLER_OUTPUT для прямого и FFTW_SCRAMBLER_INPUT для обратного БПФ.
Вычисление количества точек, обрабатываемых на каждом из процессоров:
void fftw_mpi_local_sizes(fftw_mpi_plan plan,int *local_n,int *local_start, int *local_n_after_transform, int *local_start_after_transform, int *total_local_size);
- plan — ранее созданный план
- local_n — количество точек на процессоре перед преобразованием;
- local_start — указывает с какой точки начинается нумерация на данном процессоре перед преобразованием; (таким образом, на каждом процессоре перед преобразованием размещаются точки от local_start до local_start+local_n-1)
- local_n_after_transform — количество точек на процессоре после преобразования;
- local_start_after_transform — указывает с какой точки начинается нумерация на данном процессоре после преобразования;
- total_local_size — размер рабочей области на процессоре
Функция параллельного выполнения БПФ:
fftw_mpi(fftw_mpi_plan pinv, n_fields, fftw_complex *local_data, fftw_complex *work)
Здесь n_fields — количество одинаковых векторов, которые надо преобразовать (в простейшем случае равен 1). Для выполнения многомерных преобразований в параллельной версии используются аналогичные функции с префиксом fftwnd_mpi.
Источник: parallel.ru
Функция numpy fft.fft() в Python: преобразование Фурье
Преобразование Фурье разлагает функцию на составляющие ее частоты. Частным случаем является выражение музыкального аккорда через громкость и частоту составляющих его нот.
Преобразование Фурье — это прикладная концепция в мире науки и цифровой обработки сигналов. FT (преобразование Фурье) обеспечивает представление исходного сигнала в частотной области. Например, учитывая синусоидальный сигнал во временной области, преобразование Фурье обеспечивает составляющую частоту сигнала.
С помощью преобразования Фурье периодические и непериодические сигналы могут быть преобразованы из временной в частотную область. Но сначала давайте посмотрим, как вычислять преобразования Фурье в Python.
Что такое функция numpy fft.fft() в Python?
Numpy fft.fft() — это функция, которая вычисляет одномерное дискретное преобразование Фурье в Python. Метод numpy fft.fft() вычисляет одномерное дискретное n-точечное дискретное преобразование Фурье (DFT) с помощью эффективного алгоритма быстрого преобразования Фурье (БПФ — англ. FFT) [CT].
Если вы уже установили numpy и scipy и хотите создать простое FFT набора данных, вы можете использовать функцию numpy fft.fft().
Источник: python-lab.ru