Создание спектрограммы для последовательности двухмерных видеокадров

У меня есть некоторые данные, которые состоят из последовательности видеокадров, которые представляют изменения яркости во времени относительно движущейся базовой линии. В этих видеороликах есть два типа событий, которые могут произойти - "локализованные" события, которые состоят из изменений яркости в небольших группах кластерных пикселей и загрязняющих "диффузных" событий, которые влияют на большинство пикселей в кадре:

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

Я могу найти много информации о создании спектрограмм для 1D-данных (например, аудио), но я не очень много разбираюсь в создании спектрограмм для 2D-данных. То, что я пробовал до сих пор, состоит в том, чтобы генерировать 2D-спектр мощности из преобразования Фурье кадра, затем выполнить полярное преобразование вокруг компонента постоянного тока, а затем усреднить по углам, чтобы получить 1D-спектр мощности:

Затем я применяю это для каждого кадра в моем фильме и генерирую растровый график спектральной мощности с течением времени:

Это похоже на разумный подход? Существует ли более "стандартный" подход к спектральному анализу 2D-данных?

Здесь мой код:

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

Ответ 1