为一系列 2D 电影帧生成频谱图

Posted

技术标签:

【中文标题】为一系列 2D 电影帧生成频谱图【英文标题】:Generating a spectrogram for a sequence of 2D movie frames 【发布时间】:2013-12-18 02:05:53 【问题描述】:

我有一些由一系列视频帧组成的数据,这些视频帧表示相对于移动基线的亮度随时间的变化。在这些视频中,可能会发生两种“事件”——“局部化”事件,包括一小组聚集像素的亮度变化,以及影响帧中大部分像素的污染“漫反射”事件:

我希望能够将亮度的局部变化与漫反射事件隔离开来。我计划通过减去每帧的适当低通滤波版本来做到这一点。为了设计最佳过滤器,我想知道我的帧的哪些空间频率在漫反射和局部事件期间被调制,即我想随着时间的推移生成我的电影的频谱图。

我可以找到很多关于为 1D 数据(例如音频)生成频谱图的信息,但在为 2D 数据生成频谱图方面我没有遇到太多。到目前为止,我尝试的是从帧的傅里叶变换生成 2D 功率谱,然后对 DC 分量执行极坐标变换,然后跨角度平均以获得 1D 功率谱:

然后,我将它应用到电影中的每一帧,并生成光谱功率随时间变化的光栅图:

这看起来是一种明智的做法吗?是否有更“标准”的方法来对二维数据进行光谱分析?

这是我的代码:

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

【问题讨论】:

您所做的与数字半色调世界中已知的 RAPS(径向平均功率谱)类似,因此在一般意义上它确实有意义,不确定您的视频应用程序。 .. @Jaime 知道这一点很有用 - 我会阅读 RAPS 文献 我会通过对每个帧进行阈值处理来解决这个问题,以便仅对事件进行二值化,用 ndimage 标记二值化帧并检查数字/簇大小作为上限。 @Dschoni 从那个特定的例子中并不太清楚,但也有“局部”事件同时出现在大量圆形斑点(实际上是脑细胞)中,而“扩散”事件仍然存在“模糊”,但会影响帧的空间受限部分。出于这个原因,在我的情况下,对每帧的平均像素值进行阈值处理(这是我认为你所建议的)可能不会很好地工作。我确实有一个半工作解决方案,它使用时空 ICA 来提取空间和时间上都稀疏的信号。 @ali_m:那么你如何从根本上区分“事件”和噪音?据我了解,您以某种方式(通过改变灰度值?)然后想要区分不同类型的事件。 【参考方案1】:

我相信您描述的方法通常是进行此分析的最佳方法。

但是,我确实在您的代码中发现了一个错误。如:

np.abs(f2*f2)

不是复数数组 f2 的 PSD,你需要将 f2 乘以它的复共轭而不是它本身(|f2^2| 与 |f2|^2 不同)。

相反,您应该做类似的事情

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

或者,更简洁:

np.abs(f2)**2.

2D 功率谱中的振荡是这种错误的一个明显迹象(我自己之前做过这个!)

【讨论】:

以上是关于为一系列 2D 电影帧生成频谱图的主要内容,如果未能解决你的问题,请参考以下文章

如何有效地生成和连接频谱图

FFmpeg:流式音频播放列表,标准化响度并生成频谱图和波形

如何在频谱图Python中找到峰值[重复]

在 MATLAB 中将频谱图另存为图像

关于用MATLAB设计确定信号的频谱分析和滤波

Facebook频谱图模型生成比尔·盖茨声音,性能完胜WaveNetMAESTRO