Настройка частоты синусоидальной волны

Question

Этот вопрос основан на предыдущий аналогичный вопрос.

У меня есть следующее уравнение и скорректированное (некоторые случайные данные):0,44*sin(N*2*PI/30)

Я пытаюсь использовать БПФ, чтобы получить частоту из сгенерированных данных.Однако частота оказывается близкой к частоте, но не равной ей (что делает волну немного больше, чем предполагалось).

Максимальная частота для БПФ составляет 7 Гц, однако ожидаемая частота (30/2PI) составляет 4,77 Гц.

Я включил график БПФ и нанес значения.

Код, который я использую:

[sampleFFTValues sFreq] = positiveFFT(sampledata, 1);
sampleFFTValues = abs(sampleFFTValues);
[v sFFTV]= max(sampleFFTValues)

Положительное БПФ может быть нашел здесь.По сути, он центрирует график БПФ и отсекает отрицательные сигналы.

Мой вопрос: как я могу сделать БПФ более точным, не прибегая к методу наименьших квадратов только для частоты?

Answer 1

Я не думаю, что БПФ подходит для измерения частоты с высоким разрешением (квази)периодических сигналов - см. ниже.

Каждое дискретное БПФ имеет расширение на нецелочисленных частотах элемента разрешения (то есть на любой частоте, которая не совсем соответствует одному из шагов частоты конкретного БПФ);эти «промежуточные» частоты будут размазаны/распределены вокруг ближайшего целочисленного интервала.Форма этого расширения («функция расширения») зависит от оконной функции, используемой для БПФ.Эта функция распространения – для упрощения и обобщения – либо очень узкая, но очень неровная (очень высокие пики/очень низкие впадины), либо шире, но менее неровная.Теоретически вы могли бы выполнить очень точную развертку по частоте синусоидальных волн и вычислить БПФ для каждой из них, а затем «откалибровать» форму и поведение функции, сохранив выходные данные всех БПФ вместе с частотой, которая привела к этому результату. а затем, сравнивая выходной сигнал БПФ измеряемого сигнала с ранее сохраненными результатами и находя «ближайший» результат, найдите более точную частоту.

Много усилий.

Но не делайте этого, если вам нужно измерить только частоту одного сигнала.

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

Пример:Частота дискретизации 48000 Гц, сигнал 4,77 Гц дает разрешение ~0,0005 Гц, просто измеряя длину один цикл с самым грубым подходом.(Если вы возьмете н циклов разрешение по частоте умножается на н также.)

Answer 2

Как уже упоминалось другими, вы неправильно интерпретируете частоту сигнала.Позвольте мне привести пример, чтобы прояснить несколько вещей:

Fs = 200;                        %# sampling rate
t = 0:1/Fs:1-1/Fs;               %# time vector of 1 second 
f = 6;                           %# frequency of signal
x = 0.44*sin(2*pi*f*t);          %# sine wave

N = length(x);                   %# length of signal
nfft = N;                        %# n-point DFT, by default nfft=length(x)
                                 %# (Note: it is faster if nfft is a power of 2)
X = abs(fft(x,nfft)).^2 / nfft;  %# square of the magnitude of FFT

cutOff = ceil((nfft+1)/2);       %# nyquist frequency
X = X(1:cutOff);                 %# FFT is symmetric, take first half
X(2:end -1) = 2 * X(2:end -1);   %# compensate for the energy of the other half
fr = (0:cutOff-1)*Fs/nfft;       %# frequency vector

subplot(211), plot(t, x)
title('Signal (Time Domain)')
xlabel('Time (sec)'), ylabel('Amplitude')

subplot(212), stem(fr, X)
title('Power Spectrum (Frequency Domain)')
xlabel('Frequency (Hz)'), ylabel('Power')

Теперь вы можете видеть, что пик БПФ соответствует исходной частоте сигнала 6 Гц.

[v idx] = max(X);
fr(idx)
ans = 
      6

Мы даже можем это проверить Теорема Парсеваля держит:

( sum(x.^2) - sum(X) )/nfft < 1e-6

---

Вариант 2

В качестве альтернативы мы можем использовать функции панели инструментов обработки сигналов:

%# estimate the power spectral density (PSD) using the periodogram
h = spectrum.periodogram;
hopts = psdopts(h);
set(hopts, 'Fs',Fs, 'NFFT',nfft, 'SpectrumType','onesided')

hpsd = psd(h, x, hopts);
figure, plot(hpsd)

Pxx = hpsd.Data;
fr = hpsd.Frequencies;
[v idx]= max(Pxx)
fr(idx)

avgpower(hpsd)

Обратите внимание, что эта функция использует логарифмический масштаб: plot(fr,10*log10(Pxx)) вместо plot(fr,Pxx)

Answer 3

Предполагая, что N — это время в секундах, ваша частота равна 1/30 Гц (y=A * sin( 2* PI * f * t))

Разрешение по частоте = частота дискретизации / точки БПФ

Частота выборки определяется критерием Найквиста, частота выборки (выборок в секунду) должна быть как минимум в два раза больше максимальной анализируемой частоты, например48 кГц для анализа до 24 кГц.(Для «реальных» данных хорошо иметь немного буфера).

Итак, вам может потребоваться увеличить размер вашего БПФ.

Answer 4

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

Для получения точной синусоиды при нулевом шуме забудьте о БПФ и просто решите уравнение с 3 неизвестными, что может включать всего лишь 3 или 4 точки выборки без алиасинга, алгоритмы для этого здесь и здесь.

Я перечисляю несколько других методов оценки частоты на моем Веб-страница DSP.

Answer 5

Если вы генерируете функцию, а не работаете с выборками, вы можете сгенерировать МНОГО точек и выполнить БОЛЬШОЕ БПФ, поэтому элементы разрешения по частоте будут очень маленькими для высокой точности.Но это не решит основную проблему.

Answer 6

Во-первых, поправка к вашему вопросу:(30/2PI) — это не частота.Частота вашего сигнала равна 1/30 * любой частоты дискретизации, которую вы использовали.Во-вторых, можете ли вы сказать мне, какова была длина вектора выборочных данных?Когда FFT возвращает вектор значений, значение ITH будет соответствовать f_i = i/n, где n-длина вектора, а i in [0, n-1] вы хотите, чтобы i/n ровно одинаково 1/30 для некоторого целого числа я.Другими словами, N должно равняться 30*i, т. е. N должно быть кратно 30. Итак, была ли длина использованного вами вектора кратна 30? Если нет, попробуйте сделать это, и это должно решить проблему.

Answer 7

Попробуйте оконная функция?