Вот мой ответ, как и обещал.
Первый или ваши вопросы, связанные с тем, почему вы «получаете странные результаты при увеличении частоты до 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 равным количеству выборок в сигнале.Однако использование этого формата помогает демистифицировать ваше мышление, и это становится очень важным в моем следующем примере.
** ПОБОЧНОЕ ПРИМЕЧАНИЕ: вы используете real (F) в своем примере.Если у вас нет очень веской причины извлекать только действительную часть результата БПФ, тогда гораздо чаще извлекать величину БПФ с помощью abs (F).Это эквивалент sqrt (real (F). ^ 2 + imag (F). ^ 2). **
Большую часть времени вы захотите использовать более короткую NFFT.Это может быть связано с тем, что вы, возможно, выполняете анализ в системе реального времени, или потому, что вы хотите усреднить результат многих БПФ вместе, чтобы получить представление о среднем спектре для изменяющегося во времени сигнала, или потому, что вы хотите сравнить спектрысигналы разной длительности без потери информации.Простое использование команды fft со значением NFFT <количество элементов в вашем сигнале приведет к вычислению fft из последних точек NFFT сигнала.Это немного расточительно. </p>
Следующий пример гораздо более актуален для полезного применения.Он показывает, как вы бы разбили сигнал на блоки, а затем обработали каждый блок и усреднили результат:
%//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
Количество перекрытий, которое вы выберете, будет зависеть от типа используемого вами окна. В приведенном выше примере окно Хемминга взвешивает выборки в каждом буфере в направлении нуля от центра каждого кадра. Чтобы использовать всю информацию во входном сигнале, важно использовать некоторое перекрытие. Однако, если вы просто используете простое прямоугольное окно, перекрытие становится бессмысленным, поскольку все выборки взвешиваются одинаково. Чем больше перекрытий вы используете, тем больше обработки требуется для вычисления среднего спектра.
Надеюсь, это поможет вам понять.