Как выполнить повторную выборку / ребинирование спектра?

В 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

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

Другие идеи:

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