绘制图像的 PSD 与 x/y 轴的关系答案

【问题标题】：plot the PSD of an image vs. x/y axis绘制图像的 PSD 与 x/y 轴的关系
【发布时间】：2019-01-28 21:11:43
【问题描述】：

一位研究教授让我为几个视频生成二维空间频谱密度图。我有两个问题：

如何绘制 PSD 与 x,y 轴？
我知道如何为图像生成 PSD，但不确定如何在视频上执行相同的操作。我想为视频中的每一帧获取 PSD 并取平均值，但我在 python 中实现它时遇到了困难。

下面是我的代码

curr_dir = os.getcwd()
img = cv2.imread(curr_dir+'/test.jpg',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
mag = 20*np.log(np.abs(fshift))
plt.subplot(121), plt.imshow(img,cmap='gray')
plt.subplot(122), plt.imshow(mag,cmap='gray')
plt.show()

这会生成如下内容：

我想得到这样的东西：

非常感谢任何帮助/建议！

【问题讨论】：

您是否尝试为视频中的每一帧生成 2D 空间光谱密度图，随着时间的推移制作一种 PSD 立方体？还是每帧一维 PSD，制作一个瀑布图，其中一个轴是您的 PSD，另一个轴是时间？
@Engineero 嗨，我正在尝试为视频生成空间光谱密度图（我知道这听起来很奇怪，但这就是我的要求）。我想绘制 PSD 的 x（顺风）和 y（顺风）分量与频率的关系。

标签： python numpy signal-processing

【解决方案1】：

由于您显示了两个 1d 光谱，因此您似乎正在寻找类似以下的内容。

我们在图像中读取，沿一个轴进行傅里叶变换，然后沿另一个轴对每个 bin 中的功率求和。由于输入是实值，我们使用 rfft() 所以我们不必移动频谱，我们使用 rfftreq() 来计算每个 bin 的频率。我们绘制的结果省略了 0 频率区间（对应于基线）中有时较大的信号，以便频谱的有用部分以方便的比例显示。

#!/usr/bin/python3
import cv2
import os
import math
import matplotlib
import matplotlib.pyplot as plt
import numpy as np

curr_dir = os.getcwd()
img = cv2.imread(curr_dir+'/temp.png',0)
print( img.shape )

# Fourier Transform along the first axis

# Round up the size along this axis to an even number
n = int( math.ceil(img.shape[0] / 2.) * 2 )

# We use rfft since we are processing real values
a = np.fft.rfft(img,n, axis=0)

# Sum power along the second axis
a = a.real*a.real + a.imag*a.imag
a = a.sum(axis=1)/a.shape[1]

# Generate a list of frequencies
f = np.fft.rfftfreq(n)

# Graph it
plt.plot(f[1:],a[1:], label = 'sum of amplitudes over y vs f_x')

# Fourier Transform along the second axis

# Same steps as above
n = int( math.ceil(img.shape[1] / 2.) * 2 )

a = np.fft.rfft(img,n,axis=1)

a = a.real*a.real + a.imag*a.imag
a = a.sum(axis=0)/a.shape[0]

f = np.fft.rfftfreq(n)

plt.plot(f[1:],a[1:],  label ='sum of amplitudes over x vs f_y')

plt.ylabel( 'amplitude' )
plt.xlabel( 'frequency' )
plt.yscale( 'log' )

plt.legend()

plt.savefig( 'test_rfft.png' )
#plt.show()

将此应用于您问题中发布的照片，会产生以下结果，

【讨论】：