MATLAB FFT xaxis ограничивает путаницу и смещение
Я впервые использую функцию fft и пытаюсь построить частотный спектр простой косинус-функции:
f = cos (2 * пи *300* т)
Частота дискретизации составляет 220500. Я строю одну секунду функции f.
Вот моя попытка:
time = 1;
freq = 220500;
t = 0 : 1/freq : 1 - 1/freq;
N = length(t);
df = freq/(N*time);
F = fftshift(fft(cos(2*pi*300*t))/N);
faxis = -N/2 / time : df : (N/2-1) / time;
plot(faxis, real(F));
grid on;
xlim([-500, 500]);
Почему я получаю странные результаты при увеличении частоты до 900 Гц? Эти странные результаты можно исправить, увеличив пределы оси X, скажем, с 500 Гц до 1000 Гц. Кроме того, это правильный подход? Я заметил, что многие другие люди не использовали fftshift(X) (но я думаю, что они сделали только односторонний анализ спектра).
Спасибо.
2 ответа
Вот мой ответ, как и обещал.
Первый или ваши вопросы, связанные с тем, почему вы "получаете странные результаты при увеличении частоты до 900 Гц", связаны с функцией масштабирования графика Matlab, как описано в @Castilho. Когда вы изменяете диапазон оси X, Matlab попытается помочь вам и изменить масштаб оси Y. Если пики лежат за пределами указанного диапазона, Matlab увеличит небольшие числовые ошибки, сгенерированные в процессе. Вы можете исправить это с помощью команды ylim, если это вас беспокоит.
Тем не менее, ваш второй, более открытый вопрос "это правильный подход?" требует более глубокого обсуждения. Позвольте мне рассказать вам о том, как бы я сделал более гибкое решение для достижения вашей цели построения волны косинуса.
Вы начинаете со следующего:
time = 1;
freq = 220500;
Это сразу вызывает тревогу в моей голове. Глядя на остальную часть поста, вы, кажется, интересуетесь частотами в диапазоне менее кГц. Если это так, то эта частота дискретизации является чрезмерной, поскольку предел Найквиста (sr/2) для этой частоты выше 100 кГц. Я предполагаю, что вы хотели использовать обычную частоту дискретизации звука 22050 Гц (но я могу ошибаться здесь)?
В любом случае, ваш анализ в конечном итоге работает численно. Однако вы не помогаете себе понять, как наиболее эффективно использовать БПФ для анализа в реальных ситуациях.
Позвольте мне написать, как я это сделаю. Следующий скрипт выполняет почти то же, что и ваш скрипт, но открывает некоторый потенциал, на котором мы можем строить.,
%// These are the user parameters
durT = 1;
fs = 22050;
NFFT = durT*fs;
sigFreq = 300;
%//Calculate time axis
dt = 1/fs;
tAxis = 0:dt:(durT-dt);
%//Calculate frequency axis
df = fs/NFFT;
fAxis = 0:df:(fs-df);
%//Calculate time domain signal and convert to frequency domain
x = cos( 2*pi*sigFreq*tAxis );
F = abs( fft(x, NFFT) / NFFT );
subplot(2,1,1);
plot( fAxis, 2*F )
xlim([0 2*sigFreq])
title('single sided spectrum')
subplot(2,1,2);
plot( fAxis-fs/2, fftshift(F) )
xlim([-2*sigFreq 2*sigFreq])
title('whole fft-shifted spectrum')
Вы вычисляете ось времени и вычисляете количество точек БПФ по длине оси времени. Это очень странно. Проблема этого подхода заключается в том, что разрешение по частоте БПФ изменяется при изменении длительности входного сигнала, поскольку N зависит от вашей переменной времени. Команда matlab fft будет использовать размер FFT, который соответствует размеру входного сигнала.
В моем примере я вычисляю ось частоты непосредственно из NFFT. Это несколько не имеет значения в контексте приведенного выше примера, так как я установил NFFT равным количеству выборок в сигнале. Однако использование этого формата помогает демистифицировать ваше мышление и становится очень важным в моем следующем примере.
** Боковое примечание: вы используете реальный (F) в вашем примере. Если у вас нет очень веской причины извлекать только реальную часть результата БПФ, тогда гораздо чаще извлекать величину БПФ с помощью abs(F). Это эквивалент sqrt(real(F).^2 + imag(F).^2).**
Большую часть времени вы захотите использовать более короткий NFFT. Это может быть связано с тем, что вы, возможно, выполняете анализ в системе реального времени, или потому, что вы хотите усреднить результат многих БПФ вместе, чтобы получить представление о среднем спектре для изменяющегося во времени сигнала, или потому, что вы хотите сравнить спектры сигналы разной длительности без потери информации. Простое использование команды fft со значением NFFT <количество элементов в вашем сигнале приведет к вычислению fft из последних точек NFFT сигнала. Это немного расточительно.
Следующий пример гораздо более актуален для полезного приложения. Он показывает, как вы бы разбили сигнал на блоки, а затем обработали каждый блок и усреднили результат:
%//These are the user parameters
durT = 1;
fs = 22050;
NFFT = 2048;
sigFreq = 300;
%//Calculate time axis
dt = 1/fs;
tAxis = dt:dt:(durT-dt);
%//Calculate frequency axis
df = fs/NFFT;
fAxis = 0:df:(fs-df);
%//Calculate time domain signal
x = cos( 2*pi*sigFreq*tAxis );
%//Buffer it and window
win = hamming(NFFT);%//chose window type based on your application
x = buffer(x, NFFT, NFFT/2); %// 50% overlap between frames in this instance
x = x(:, 2:end-1); %//optional step to remove zero padded frames
x = ( x' * diag(win) )'; %//efficiently window each frame using matrix algebra
%// Calculate mean FFT
F = abs( fft(x, NFFT) / sum(win) );
F = mean(F,2);
subplot(2,1,1);
plot( fAxis, 2*F )
xlim([0 2*sigFreq])
title('single sided spectrum')
subplot(2,1,2);
plot( fAxis-fs/2, fftshift(F) )
xlim([-2*sigFreq 2*sigFreq])
title('whole fft-shifted spectrum')
Я использую окно Хэмминга в приведенном выше примере. Выбранное вами окно должно соответствовать приложению http://en.wikipedia.org/wiki/Window_function
Выбранная вами величина перекрытия будет в некоторой степени зависеть от типа используемого вами окна. В приведенном выше примере окно Хэмминга взвешивает выборки в каждом буфере в направлении нуля от центра каждого кадра. Чтобы использовать всю информацию во входном сигнале, важно использовать некоторое перекрытие. Однако, если вы просто используете простое прямоугольное окно, перекрытие становится бессмысленным, так как все выборки взвешиваются одинаково. Чем больше перекрытий вы используете, тем больше обработки требуется для вычисления среднего спектра.
Надеюсь, это поможет вам понять.
Ваш результат совершенно прав. Ваш расчет частотной оси тоже правильный. Проблема лежит на шкале оси Y. Когда вы используете функцию xlims, matlab автоматически пересчитывает шкалу y, чтобы вы могли видеть "значимые" данные. Когда пики косинуса лежат за пределами выбранного вами предела (при f>500 Гц), пики не отображаются, поэтому масштаб рассчитывается на основе некоторого очень маленького шума (здесь, на моем компьютере, с matlab 2011a, масштаб y был 10-16).
Изменение предела действительно является правильным подходом, потому что если вы не измените его, вы не сможете увидеть пики на частотном спектре.
Однако я заметил одну вещь. Есть ли причина для вас, чтобы представить реальную часть преобразования? Обычно это abs(F) это заговор, а не реальная часть.
редактировать: На самом деле, ваша ось частоты верна только потому, что df, в данном случае, равен 1. Факсимильная линия верна, а вычисление df - нет.
БПФ рассчитывает N точек от -Fs/2 до Fs/2. Таким образом, N точек в диапазоне Fs дает df Fs/N. Как N/ время = Fs => время = N/Fs. Подставляя это в выражение df, вы использовали: your_df = Fs/N*(N/Fs) = (Fs/N)^2. Поскольку Fs/N = 1, конечный результат оказался верным:P