Преобразование Фурье с помощью Matlab

Ответ 1

TL; DR (резюме)

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

Взять домой сообщения являются:

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

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

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

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

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

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

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

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

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

Код

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

Ответ 2

Я нашел этот онлайн пример fft (https://habr.com/post/112068/) весьма полезным. Проверь это:

%% Parameters
Tm=5;% Length of signal (s)
Fd=512;% Sampling frequency (Hz)
Ak=0.5;% Constant component (Unit)
A1=1;% The amplitude of the first sinusoid (Unit)
A2=0.7;% Amplitude of the second sinusoid (Unit)
F1=13;% Frequency of the first sinusoid (Hz)
F2=42;% Frequency of the second sinusoid (Hz)
Phi1=0;% Initial phase of the first sinusoid (Degrees)
Phi2=37;% The initial phase of the second sinusoid (Degrees)
An=3*A1;% Noise Dispersion (Unit)
FftL=1024;% Number of Fourier Spectrum Lines
%% Generating work arrays
T=0:1/Fd:Tm;% Time Arrays
Noise=An*randn(1,length(T));% An array of random noise with a length equal to the array of time
Signal=Ak+A1*sind((F1*360).*T+Phi1)+A2*sind((F2*360).*T+Phi2);% Signal array (a mixture of 2x sinusoids and a constant component)
%% Spectral representation of the signal
FftS=abs(fft(Signal,FftL));% The amplitudes of the Fourier transform of the signal
FftS=2*FftS./FftL;% Spectrum normalization by amplitude
FftS(1)=FftS(1)/2;% The normalization of the constant component in the spectrum
FftSh=abs(fft(Signal+Noise,FftL));% The amplitudes of the Fourier transform of the signal + noise mixture
FftSh=2*FftSh./FftL;% Spectrum normalization by amplitude
FftSh(1)=FftSh(1)/2;% The normalization of the constant component in the spectrum
%% Plotting
subplot(2,1,1);
plot(T,Signal);
title('Signal');
xlabel('Time (s)');
ylabel('Amplitude (Unit)');
subplot(2,1,2);
plot(T,Signal+Noise);
title('Signal+Noise');
xlabel('Time (s)');
ylabel('Amplitude (Unit)');

F=0:Fd/FftL:Fd/2-1/FftL;% The frequency array of the calculated Fourier spectrum
figure();
subplot(2,1,1);
plot(F,FftS(1:length(F)));%  Plotting of the spectrum of the Fourier signal
title('Signal spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude (Unit)');
subplot(2,1,2);
plot(F,FftSh(1:length(F)));% Plotting of the Fourier signal spectrum
title('Signal spectrum');
xlabel('Frequency (Hz)');
ylabel('Amplitude (Unit)');