image-processing - señales - espectrograma pdf

Generación de un espectrograma para una secuencia de cuadros de película 2D

Question

Tengo algunos datos que consisten en una secuencia de cuadros de video que representan cambios en la luminancia a lo largo del tiempo en relación con una línea base en movimiento. En estos videos hay dos tipos de ''eventos'' que pueden ocurrir: eventos ''localizados'', que consisten en cambios de luminancia en pequeños grupos de píxeles agrupados y eventos contaminantes ''difusos'', que afectan a la mayoría de los píxeles en el marco:

Me gustaría poder aislar los cambios locales en la luminancia de los eventos difusos. Estoy planeando hacerlo restando una versión filtrada de paso bajo apropiada de cada fotograma. Para diseñar un filtro óptimo, me gustaría saber qué frecuencias espaciales de mis cuadros se modulan durante los eventos difusos y locales, es decir, me gustaría generar un espectrograma de mi película a lo largo del tiempo.

Puedo encontrar mucha información sobre la generación de espectrogramas para datos 1D (por ejemplo, audio), pero no he encontrado mucho en la generación de espectrogramas para datos 2D. Lo que he intentado hasta ahora es generar un espectro de potencia 2D a partir de la transformada de Fourier del cuadro, luego realizar una transformación polar sobre el componente CC y promediar en todos los ángulos para obtener un espectro de potencia 1D:

A continuación, aplico esto a cada fotograma de mi película y genero un gráfico de trama de potencia espectral a lo largo del tiempo:

¿Esto parece un enfoque sensato para tomar? ¿Hay un enfoque más "estándar" para hacer un análisis espectral en datos 2D?

Aquí está mi código:

import numpy as np
# from pyfftw.interfaces.scipy_fftpack import fft2, fftshift, fftfreq
from scipy.fftpack import fft2, fftshift, fftfreq
from matplotlib import pyplot as pp
from matplotlib.colors import LogNorm
from scipy.signal import windows
from scipy.ndimage.interpolation import map_coordinates

def compute_2d_psd(img, doplot=True, winfun=windows.hamming, winfunargs={}):

    nr, nc = img.shape
    win = make2DWindow((nr, nc), winfun, **winfunargs)

    f2 = fftshift(fft2(img*win))
    psd = np.abs(f2*f2)
    pol_psd = polar_transform(psd, centre=(nr//2, nc//2))

    mpow = np.nanmean(pol_psd, 0)
    stdpow = np.nanstd(pol_psd, 0)

    freq_r = fftshift(fftfreq(nr))
    freq_c = fftshift(fftfreq(nc))
    pos_freq = np.linspace(0, np.hypot(freq_r[-1], freq_c[-1]), 
        pol_psd.shape[1])

    if doplot:
        fig,ax = pp.subplots(2,2)

        im0 = ax[0,0].imshow(img*win, cmap=pp.cm.gray)
        ax[0,0].set_axis_off()
        ax[0,0].set_title(''Windowed image'')

        lnorm = LogNorm(vmin=psd.min(), vmax=psd.max())
        ax[0,1].set_axis_bgcolor(''k'')
        im1 = ax[0,1].imshow(psd, extent=(freq_c[0], freq_c[-1], 
            freq_r[0], freq_r[-1]), aspect=''auto'', 
            cmap=pp.cm.hot, norm=lnorm)
        # cb1 = pp.colorbar(im1, ax=ax[0,1], use_gridspec=True)
        # cb1.set_label(''Power (A.U.)'')
        ax[0,1].set_title(''2D power spectrum'')

        ax[1,0].set_axis_bgcolor(''k'')
        im2 = ax[1,0].imshow(pol_psd, cmap=pp.cm.hot, norm=lnorm, 
            extent=(pos_freq[0],pos_freq[-1],0,360), 
            aspect=''auto'')
        ax[1,0].set_ylabel(''Angle (deg)'')
        ax[1,0].set_xlabel(''Frequency (cycles/px)'')
        # cb2 = pp.colorbar(im2, ax=(ax[0,1],ax[1,1]), use_gridspec=True)
        # cb2.set_label(''Power (A.U.)'')
        ax[1,0].set_title(''Polar-transformed power spectrum'')

        ax[1,1].hold(True)
        # ax[1,1].fill_between(pos_freq, mpow - stdpow, mpow + stdpow, 
        #   color=''r'', alpha=0.3)
        ax[1,1].axvline(0, c=''k'', ls=''--'', alpha=0.3)
        ax[1,1].plot(pos_freq, mpow, lw=3, c=''r'')
        ax[1,1].set_xlabel(''Frequency (cycles/px)'')
        ax[1,1].set_ylabel(''Power (A.U.)'')
        ax[1,1].set_yscale(''log'')
        ax[1,1].set_xlim(-0.05, None)
        ax[1,1].set_title(''1D power spectrum'')

        fig.tight_layout()

    return mpow, stdpow, pos_freq

def make2DWindow(shape,winfunc,*args,**kwargs):
    assert callable(winfunc)
    r,c = shape
    rvec = winfunc(r,*args,**kwargs)
    cvec = winfunc(c,*args,**kwargs)
    return np.outer(rvec,cvec)

def polar_transform(image, centre=(0,0), n_angles=None, n_radii=None):
    """
    Polar transformation of an image about the specified centre coordinate
    """
    shape = image.shape
    if n_angles is None:
        n_angles = shape[0]
    if n_radii is None:
        n_radii = shape[1]
    theta = -np.linspace(0, 2*np.pi, n_angles, endpoint=False).reshape(-1,1)
    d = np.hypot(shape[0]-centre[0], shape[1]-centre[1])
    radius = np.linspace(0, d, n_radii).reshape(1,-1)
    x = radius * np.sin(theta) + centre[0]
    y = radius * np.cos(theta) + centre[1]

    # nb: map_coordinates can give crazy negative values using higher order
    # interpolation, which introduce nans when you take the log later on
    output = map_coordinates(image, [x, y], order=1, cval=np.nan, 
        prefilter=True)
    return output

Answer 1

Creo que el enfoque que describes es en general la mejor forma de hacer este análisis.

Sin embargo, detecté un error en tu código. como:

np.abs(f2*f2)

no es la PSD de la matriz compleja f2, necesita multiplicar f2 por su complejo conjugado en lugar de por sí mismo (| f2 ^ 2 | no es lo mismo que | f2 | ^ 2).

En cambio, deberías hacer algo como

(f2*np.conjugate(f2)).astype(float)

O, más limpiamente:

np.abs(f2)**2.

Las oscilaciones en el espectro de potencia 2D son un signo revelador de este tipo de error (¡he hecho esto antes que yo!)