二维功率谱绘制：频率向量生成与波数范围疑问

Great question—working with 2D FFT frequency axes can feel counterintuitive when you’re coming from 1D signal processing, since images don’t come with an explicit "sampling rate" label like audio or time-series data. Let’s unpack this step by step.

First, let’s clarify what the frequency axes represent for an image. When you run fft2 and fftshift, you’re shifting the zero-frequency component to the center of the spectrum. To map each pixel in the spectrum to a real-world frequency, you need to account for:

• The image dimensions [M, N] (rows x columns)
• The spatial sampling rate (more on this next)

Your initial attempt at fx = [fs/2*linspace(0,1,N/2+1)] misses the negative frequency half (since fftshift centers the spectrum, we need frequencies ranging from -fs/2 to fs/2). Here’s the right way to generate full 2D frequency vectors in MATLAB:

For a [M, N] image: #

[M, N] = size(img);
fs_x = 1; % We’ll explain fs_x/fs_y below
fs_y = 1;

% Generate horizontal (x/column) frequency vector
fx = linspace(-fs_x/2, fs_x/2 - fs_x/N, N);
% Generate vertical (y/row) frequency vector
fy = linspace(-fs_y/2, fs_y/2 - fs_y/M, M);

Alternatively, you can use index-based calculation for clarity:

fx = (-N/2 : N/2 - 1) / N * fs_x;
fy = (-M/2 : M/2 - 1) / M * fs_y;

Both approaches give you a symmetric frequency axis centered at 0, matching the shifted FFT output from fftshift.

The spatial sampling rate fs (often written as f_sx for horizontal, f_sy for vertical) is the reciprocal of the physical distance between adjacent pixels. Here’s how to choose it:

• No physical scale info: If you’re working with a standard digital image (like cameraman.tif) and don’t know the real-world pixel spacing, default to fs_x = fs_y = 1. This gives you frequencies in cycles per pixel, which is the standard for general image processing tasks.
• With physical scale info: If you know, for example, that each pixel corresponds to 0.005 mm (5 microns), then fs_x = 1 / 0.005 = 200 cycles/mm. This converts your frequency axes to real-world physical units.

The k₁/k₂ axes in the paper are spatial wave numbers, which are closely related to frequency. Here’s the breakdown:

• Wave number k is typically defined as k = 2πf (angular wave number, units of radians per unit distance) or sometimes k = f (cycles per unit distance, depending on field conventions).
• A range of -200 to 200 suggests the paper’s Nyquist frequency is 200. This could come from:
  • A physical sampling rate of fs = 400 cycles/unit (so Nyquist is fs/2 = 200), if they’re using k = f.
  • Another common scenario: They normalized the frequency axes by multiplying by half the image dimension. For example, if their image is 400x400 pixels and fs=1, the raw frequency range is -0.5 to 0.5 cycles/pixel. Multiplying by 400 gives -200 to 200.

Here’s a complete snippet to compute and plot the 2D spectrum with correct frequency axes:

img = imread('cameraman.tif');
[M, N] = size(img);

% Use default sampling rate (cycles per pixel)
fs_x = 1;
fs_y = 1;

% Generate frequency vectors
fx = linspace(-fs_x/2, fs_x/2 - fs_x/N, N);
fy = linspace(-fs_y/2, fs_y/2 - fs_y/M, M);

% Compute FFT and magnitude (log scale for better visibility)
F = fftshift(fft2(double(img)));
F_mag = log10(abs(F) + 1); % Add 1 to avoid log(0)

% Plot the spectrum
figure;
imagesc(fx, fy, F_mag);
xlabel('Horizontal Frequency (cycles/pixel)');
ylabel('Vertical Frequency (cycles/pixel)');
colorbar;
title('2D Fourier Spectrum of Cameraman Image');
axis xy; % Fix y-axis to start at bottom (matches image coordinates)

内容的提问来源于stack exchange，提问作者Gregor Isack