Как сделать преобразование фурье в матлаб

Обновлено: 06.07.2024

Очень сложно найти онлайн-примеры использования fft с Matlab, которые правильно нормализуют значения амплитуды/мощности. Это имеет решающее значение, если вы собираетесь сравнивать эти значения по разным сигналам различной длины. Это чаще всего проблема для вещественных входов, потому что в этом случае обычно извлекается односторонний спектр, и поэтому изменения величины должны применяться вручную при вычислении значений амплитуды или мощности. Вы можете найти gist на GitHub здесь (сообщите мне об ошибках).

    НЕ нормализуйте коэффициенты DFT напрямую (например, не пишите Y = fft(x)/L );
    Если вы используете преобразование n-точек (например, fft(x,nfft) ), то нормализатор nfft , а не numel(x) ;
    Если вы извлекаете односторонний спектр, вам нужно отрегулировать значения амплитуды/мощности, в зависимости от которых соответствуют коэффициенты DFT сопряженных пар,
    Если вы извлекаете односторонний спектр, вы должны вычислить амплитуду и мощность отдельно (т.е. не вычислять мощность из амплитуд).

Вопрос

Мой вопрос похож на, но более общий, чем этот пост, и я думаю, что есть ошибка в отношении нормализации, с новейшей версией Matlab (2015) ) так или иначе. Я помедлил опубликовать это на CodeReview SE, если вы считаете, что это будет более подходящим, сообщите мне в комментариях.

Я хотел бы проверить следующий код преобразования Фурье с использованием Matlab fft , потому что я нашел конфликтующие источники информации в Интернете, в том числе и в самой помощи Matlab, и у меня есть не смог проверить теорему Парсеваля с некоторыми такими "рецептами" (в том числе с answers от команды MathWorks, см. ниже), особенно те извлечение односторонних спектров для реальных входов.

Например, частое удвоение амплитуды, обычно найденное в режиме онлайн для учета симметричных спектров вещественных сигналов при извлечении положительных частот, только кажется неправильным (теорема Парсеваля не выполняется), и вместо этого представляется необходимым использовать квадрат -root двух коэффициентов в Matlab (я не знаю почему). Некоторые люди также, похоже, нормализуют коэффициенты DFT напрямую, как Y = fft(x)/L , но я думаю, что это запутывает и должно быть обескуражено; амплитуда определена как модуль сложного коэффициента DFT, деленный на длину сигнала, сами коэффициенты не должны быть разделены. После проверки я намерен опубликовать этот код как сущность GitHub.

Пример использования

Неверный пример 1

Проблема и решение:

Приходит из нормализации в строке Y = fft(y,NFFT)/L .
Это должно быть вместо:

Неверный пример 2

Из команды MathWorks собственный пояснительный отчет:

Проблема и решение:

Нормализация снова. Замените Y = fft(y,NFFT)/L; на Y = fft(y,NFFT) и предположительно MX=2*abs(Y); на MX=2*abs(Y)/NFFT; . Но здесь возникает проблема удвоения амплитуды; поправочный коэффициент представляется sqrt(2) , а не 2 .

Неверный пример 3

Найдено как answer на MatlabCentral:

Проблема и решение:

Как и в первом примере, проблема нормализации. Вместо этого напишите:

Думаю, я выяснил две основные проблемы; оба имеют отношение к тому, как нормализуются коэффициенты. Позвольте мне повторить еще раз, я не мог бы препятствовать вам достаточно НЕ, чтобы напрямую нормализовать сложные коэффициенты DFT, возвращаемые функцией fft .

В дальнейшем допустим следующую номенклатуру:

Амплитудные, силовые и односторонние спектры

Как определено и объяснено в Википедии:

    Коэффициенты DFT являются сложными и не нормализованы, а формула для обратного ДПФ несет фактор 1/N перед суммой. Это естественно в некотором смысле, так как перемещение во времени к частоте может рассматриваться как проекция на базис (ортогональных) волн с разными частотами, тогда как перемещение в направлении "время от времени" можно рассматривать как взвешенную суперпозицию волн.
    В этой суперпозиции общая величина должна быть равна исходной временной точке (т.е. инверсии), а амплитуда каждой волны в этом средневзвешенном значении представляет собой просто величину соответствующий коэффициент DFT, деленный на количество волн |Xk| / N . Точно так же мощность каждой волны |Xk|^2 / N . Matlab также использует эту нормировку (ну, FFTW делает).
    Для вещественных входов коэффициенты DFT являются сопряженными парами для соответствующих положительных/отрицательных частот, кроме компонент постоянного тока (средний член, частота 0) и частоту Найквиста, когда количество точек времени четное. На практике большинство приложений избегают этой избыточности, извлекая коэффициенты DFT только для положительных частот. Это приводит к усложнениям дискретных значений амплитуды и мощности.
    Амплитуды, соответствующие парным коэффициентам DFT (все, кроме первого и частота Найквиста, когда он существует), могут быть просто удвоены, а отрицательная частота затем отбрасывается. То же самое для мощности.
    Аналогично для мощности, но обратите внимание, что в этом случае неверно для вычисления дискретных значений мощности с использованием скорректированных значений амплитуды. Действительно, N * amp_adjusted[k]^2 = N * (2*|Xk|/N)^2 не равно 2*|Xk|^2 / N (здесь квадратный корень из двух пришел в OP). Поэтому необходимо вычислить независимо значения амплитуды и мощности из коэффициентов DFT (еще одна веская причина - не масштабировать их).

Преобразование N-точек

Многие из примеров онлайн используют явное преобразование N-точек: Y = fft(x,NFFT) , где NFFT обычно имеет мощность 2, что делает вычисления более эффективными с FFTW.

Эффективная разница в этом случае (при условии, что NFFT >= N ) заключается в том, что x заполняется 0 на своем конце, пока не достигнет длины точек времени NFFT. Это означает, что число частот в разложении изменяется, и поэтому нормализация должна выполняться относительно NFFT волновых компонентов, а не исходных N временных точек.

Следовательно, почти все примеры, найденные в Интернете, неверны в том, как они нормализуют коэффициенты. Это не должно быть Y = fft(x,NFFT)/N , но Y = fft(x,NFFT)/NFFT - еще одна веская причина потерять привычку нормализовать комплексные коэффициенты.

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

Поэтому код в OP неверен, и вместо этого возникает необходимость вывести как амплитуду, так и мощность, так как нет общего коэффициента нормализации, который мог бы учитывать сложные и реальные случаи с четным или нечетным числом раз -точки. Здесь вы можете найти Gist .

Читайте также: