如何通过径向平均从二维噪声功率谱计算一维功率谱答案

【问题标题】：How to calculate 1D power spectrum from 2D noise power spectrum by radial averaging如何通过径向平均从二维噪声功率谱计算一维功率谱
【发布时间】：2015-05-24 13:50:46
【问题描述】：

伙计们，我正在尝试从图像的 2D FFT 计算 1D 功率谱。我用水平平均来做到这一点，但通过查看图表，我觉得没有意义。您能否建议如何对 2D 数据集进行径向平均以达到噪声功率谱的 1D 表示。谢谢我会感谢你的帮助。

这是我的代码 $

 fid = fopen('C:\Users\3772khobrap\Desktop\project  related\NPS_cal_data_UB\100000006.raw','r');
img = fread(fid,[512 512],'uint16');
roi = zeros(64);
avg = zeros(64);
Ux= 0.0075;% Pixel size
Uy = 0.0075;% Pixel size
%% This block of code is subdividing imaage into smaller ROI and averaging purpose
for r = 1:8 
    r_shift = (r-1)*64;
    for c = 1:8 
        c_shift = (c-1)*64;
        for i = 1:64 
            for j = 1:64
                p = img(i+r_shift,j+c_shift);
                roi(i,j) = p;
            end
        end
        avg = avg+roi;
    end
end
avg = avg./64;
%%Actual process of NPS calculation
scale = (Ux*Uy)./(64*64);%Scaling fator as per NPS calculation formula
f_x = 1/(2*Ux);%Nyquiest frequecy along x direction
f_y = 1/(2*Ux);%Nyquiest frequecy along y direction
FFT_2d = (fftshift(fft2(avg))).^2;% Power spectrum calculation
NPS = abs(FFT_2d).*scale; %% 2D NPS 
f = linspace(-f_x,f_y,64);% x-axis  vector for 1D NPS 
X_img = linspace(-f_x,f_x,512);% X axis of NPS image
Y_img = linspace(-f_x,f_x,512);% Y axis of NPS image

figure(1)
subplot(2,2,1)
imagesc(X_img,Y_img,img)
colormap gray
xlabel('X [cm]'); ylabel('Y [cm]')
title('noise image')
subplot(2,2,2)
imagesc(f,f,log(NPS))
colormap gray
xlabel('frequency [cm^{-1}]'); ylabel('frequency [cm^{-1}]');
title('2D NPS')
subplot(2,2,3)
plot(f_p,NPS(:,32))
xlabel('frequency [cm^{-2}]'); ylabel('NPS [cm^{-2}]')
title('1D NPS from central slice')
subplot(2,2,4)
plot(f_p,mean(NPS,2))
xlabel('frequency [cm^{-2}]'); ylabel('NPS [cm^2]')
title('1D NPS along X direction')

【问题讨论】：

我不明白。您应该如何从 2D 信号中获得 1D 功率谱？水平平均没有意义，因为每行的强度分布可能有很大不同，因此平均无法充分捕捉行为。
google 搜索返回了this 和this。这就是你要找的吗？
@rayreng，我想这个想法是不想区分只是方向不同但具有相同空间频率或波数的频率分量，所以频谱超过 |k| = sqrt(k1^2 + k2^2).
感谢您的回复。@rayryeng 和 A.Donda。假设对象对称的傅里叶变换，我们可以在径向网格上进行插值，并且对于每个半径，我们可以计算命中的数量。所有半径的平均值将为我们提供一维功率谱。

标签： matlab image-processing fft noise

【解决方案1】：

你可以编写这样的函数：

function profile = radialAverage(IMG, cx, cy, w)
    % computes the radial average of the image IMG around the cx,cy point
    % w is the vector of radii starting from zero
    [a,b] = size(IMG);
    [X, Y] = meshgrid( (1:a)-cx, (1:b)-cy);
    R = sqrt(X.^2 + Y.^2);
    profile = [];
    for i = w % radius of the circle
        mask = (i-1<R & R<i+1); % smooth 1 px around the radius
        values = (1-abs(R(mask)-i)) .* double(IMG(mask)); % smooth based on distance to ring
        % values = IMG(mask); % without smooth
        profile(end+1) = mean( values(:) );
    end
end

【讨论】：