Как сделать преобразование фурье в матлаб
Очень сложно найти онлайн-примеры использования 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 .
Читайте также: