¿Cómo remuestrear/rebinar un espectro?

Question

En Matlab, frecuentemente calculo los espectros de potencia usando el método de Welch (pwelch), que luego despliego en un diagrama de log-log. Las frecuencias estimadas por pwelch están igualmente espaciadas, sin embargo, los puntos espaciados logarítmicamente serían más apropiados para el gráfico log-log. En particular, al guardar la trama en un archivo PDF, esto resulta en un gran tamaño de archivo debido al exceso de puntos en alta frecuencia.

¿Cuál es un esquema efectivo para volver a muestrear (rebin) el espectro, desde frecuencias espaciadas linealmente a frecuencias espaciadas de registro? O, ¿qué es una forma de incluir espectros de alta resolución en archivos PDF sin generar tamaños de archivos excesivamente grandes?

Lo obvio que hacer es simplemente usar 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)); 

Sin embargo, esto no es deseable, ya que no conserva la energía en el espectro. Por ejemplo, si hay una gran línea espectral entre dos de los nuevos compartimientos de frecuencia, simplemente se excluirá del espectro muestreado de registro resultante.

Para solucionar este problema, nosotros en cambio podemos interpolar la integral del espectro de potencia:

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

Esto es lindo y sobre todo funciona, pero los centros de basura no son del todo bien, y no lo hace maneja de forma inteligente la región de baja frecuencia, donde la cuadrícula de frecuencia puede ser más finamente muestreada.

Otras ideas:

• escribir una función que determina la nueva frecuencia de hurgar en la basura y luego se utiliza para hacer el accumarray rebinning.
• Aplicar un filtro suavizante al espectro antes de realizar la interpolación. Problema: el tamaño del núcleo suavizado debería ser adaptable al suavizado logarítmico deseado.
• La función pwelch acepta un argumento vectorial de frecuencia f, en cuyo caso calcula el PSD a las frecuencias deseadas utilizando el algoritmo Goetzel. Tal vez solo llamar a pwelch con un vector de frecuencia espaciado en el registro sería adecuado. (¿Es esto más o menos eficiente?)
• Para el problema del tamaño del archivo PDF: incluya una imagen de mapa de bits del espectro (parece kludgy - ¡Quiero buenos gráficos vectoriales!);
• o quizás mostrar región (intervalo de confianza/polígono) en lugar de simplemente una línea segmentada para indicar el espectro.

Answer 1

Solución encontrado: https://dsp.stackexchange.com/a/2098/64

Brevemente, una solución a este problema es llevar a cabo el método de Welch con una longitud dependiente de la frecuencia de transformación. El enlace de arriba es una respuesta dsp.SE que contiene una cita en papel y una implementación de muestra. Una desventaja de esta técnica es que no puede usar la FFT, pero debido a que la cantidad de puntos DFT que se calcula se reduce considerablemente, este no es un problema grave.

Answer 2

Dejo que haga el trabajo por mí y le da las frecuencias desde el principio. El doc indica que las frecuencias que especifique se redondearán a la papelera DFT más cercana. Eso no debería ser un problema ya que está usando los resultados para trazar. Si le preocupa el tiempo de ejecución, simplemente lo probaría y mediría el tiempo.

Si quiere recargarlo usted mismo, creo que es mejor que solo escriba su propia función para hacer la integración de cada uno de sus nuevos bins. Si desea hacer su vida más fácil, puede hacer lo que hace y asegurarse de que sus compartimientos de registro compartan los límites con los lineales.

Answer 3

Si desea volver a muestrear una FFT a una tasa variable (logarítmicamente), entonces el kernel de filtro de paso bajo o suavizado necesitará ser de ancho variable para evitar el alias (pérdida de puntos de muestra). Simplemente use un kernel de interpolación de sincronización de ancho diferente para cada punto de trazado (ancho de sincronización aproximadamente el recíproco de la tasa de muestreo local).