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

В Matlab я часто вычисляю спектры мощности, используя метод Уэлча (pwelch), который я затем отображаю на графике лог-лога. Частоты, оцениваемые по pwelch, расположены на равном расстоянии друг от друга, однако логарифмически расположенные точки были бы более подходящими для логарифмического графика. В частности, при сохранении графика в 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: включите растровое изображение спектр (кажется, клуджи-я хочу хорошую векторную графику!);
  • или, возможно, отображать область (полигон/доверительный интервал) вместо простой сегментированной линии для обозначения спектра.
3 7
matlab fft spectrum

3 ответа:

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

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

2

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

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

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

0