Как пересчитать / пересобрать спектр?

Question

В Matlab я часто вычисляю спектры мощности, используя метод Уэлча (pwelch), который я затем отображаю на графике log-log. Частоты, оцененные как pwelch, расположены на одинаковом расстоянии, однако логарифмически разнесенные точки были бы более подходящими для графика log-log. В частности, при сохранении графика в файл PDF это приводит к огромному размеру файла из-за избытка точек на высокой частоте.

Какова эффективная схема для повторной выборки (ребинга) спектра, от линейно разнесенных частот до логарифмически разнесенных частот? Или как можно включить спектры высокого разрешения в файлы PDF без чрезмерного создания большие размеры файлов?

Очевидное, что нужно сделать, это просто использовать interp1:

  rate = 16384; %# sample rate (samples/sec)  
  nfft = 16384; %# number of points in the fft

  [Pxx, f] =  pwelch(detrend(data), hanning(nfft), nfft/2, nfft, rate);

  f2 = logspace(log10(f(2)), log10(f(end)), 300);
  Pxx2 = interp1(f, Pxx, f2);

  loglog(f2, sqrt(Pxx2)); 

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

Чтобы исправить это, мы можем вместо этого интерполировать интеграл спектра мощности:

  df = f(2) - f(1);
  intPxx = cumsum(Pxx) * df;                     % integrate
  intPxx2 = interp1(f, intPxx, f2);              % interpolate
  Pxx2 = diff([0 intPxx2]) ./ diff([0 F]);       % difference

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

Другие идеи:

• написать функцию, которая определяет новое частотное преобразование, а затем использует accumarray для выполнения повторного объединения.
• Примените сглаживающий фильтр к спектру перед выполнением интерполяции. Проблема: размер ядра сглаживания должен быть адаптивным к желаемому логарифмическому сглаживанию.
• Функция pwelch принимает аргумент частотного вектора f, и в этом случае она вычисляет PSD на желаемых частотах, используя алгоритм Гетцеля. Возможно, достаточно было бы просто вызвать pwelch с частотным вектором с логарифмическим интервалом. (Это более или менее эффективно?)
• Для проблемы размера файла PDF: включите растровое изображение спектра (кажется грязным - я хочу хорошую векторную графику!);
• или, возможно, отобразить регион (полигон / доверительный интервал) вместо простой сегментированной линии для обозначения спектра.

Answer 1

Я бы позволил ему сделать работу за меня и дать ему частоты с самого начала.В документе указано, что указанные вами частоты будут округлены до ближайшей корзины DFT.Это не должно быть проблемой, так как вы используете результаты для построения графика.Если вы беспокоитесь о времени выполнения, я бы просто попробовал его и рассчитал время.

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

Answer 2

Если вы хотите пересчитать FFT с переменной скоростью (логарифмически), то ядро ​​сглаживающего фильтра или фильтра нижних частот также должно быть переменной ширины, чтобы избежать наложения (потери точек выборки). Просто используйте различную ширину ядра интерполяции синхронизации для каждой точки графика (ширина синхронизации приблизительно равна обратной частоте локальной выборки).

Answer 3

Решение найдено: https://dsp.stackexchange.com/a/2098/64

Вкратце, одним из решений этой проблемы является выполнение метода Уэлча с зависимой от частоты длиной преобразования. Приведенная выше ссылка относится к ответу dsp.SE, содержащему цитирование на бумаге и пример реализации. Недостаток этого метода заключается в том, что вы не можете использовать БПФ, но поскольку число вычисляемых точек ДПФ значительно сокращено, это не является серьезной проблемой.