Why do these PSD differ?

Question

I have a circular 2D convolution (the underlying function $f$ only depends on radius) that I apply to a 2D gaussian white noise map $M$.

I am then comparing :

• the PSD (Power Spectral Density) along one axis of the 2D FFT of $f$, meaning :

a) computing the 2D FFT of $f(r)$ sampled over a significant 2D domain

b) extracting a 1D line out of this 2D FFT

• the PSD (Power Spectral Density) over a line of $f*M$, meaning :

a) computing the convoluted 2D map $f*M$

b) computing the PSD over an axis of this resulting map

And the resulting PSDs don't match. Yet it is my understanding they should.

Here is the resulting plot (Matlab source provided below) :

• First, is this normal ?
• Second, can you explain what is happening here and potentially suggest corrections ?

rng(0)

Te=0.05; % sampling time
ctn_noise_psd=1; % Bilateral PSD of requested continuous noise
sbd=sqrt(ctn_noise_psd/(2*Te)); % std of corresponding discrete noise

n1d=2000;
noise2D=randn(n1d,n1d*2)*sbd;

% Convolution to apply
% --------------------

% Continuous
tau=1;
fconv=@(r) exp(-(r/tau).^2/2);

% Numerical matrix

% we build the matrix of ranges wrt to its center
rmax=5*tau;ir=ceil(0.5*rmax/Te);rmax=2*Te*ir;
X=ones(2*ir+1,1)*linspace(0,1,2*ir+1);Y=X'; % adim : de 0 à 1
Rquart=rmax*sqrt(X.^2+Y.^2); % de 0 à rmax
R=[Rquart(end:-1:2,end:-1:2) Rquart(end:-1:2,:);Rquart(:,end:-1:2) Rquart]; % rmax around 0

% and the convolution matrix is
mconv=fconv(R);

% 2D FFT of this convolution matrix
ltf=2^(1+nextpow2(size(mconv,1)));
freqs_fft=(0:ltf/2)/(ltf*Te);
fft2D=fft2(mconv,ltf,ltf); % fft 2d

% "Cut" of the FFT along one direction
fftcut=fft2D(1:ltf/2+1,1:ltf/2+1); % 1D cut from fft2D
dsp_fft=Te^4*abs(fftcut(1,:)).^2; % corresponding PSD


% Application of the convolution
% ------------------------------
irp=ceil(5*tau/Te);
%residu=conv2fft(noise2D,mconv,'same');
residu=conv2(noise2D,mconv,'same');
residu=residu(irp+1:end-irp,irp+1:end-irp);
[nl,nc]=size(residu);

% Frequency analysis of the result
% 
[Pxx,w]=pwelch(residu(1,:));
ratio=2*pi*Te;
freqs_ech=w'/ratio;
dsp_echs=zeros(nl,length(freqs_ech));
dsp_echs(1,:)=Pxx'*ratio;
for kl=2:nl
    dsp_echs(kl,:)=pwelch(residu(kl,:))'*ratio;
end
dsp_ech=mean(dsp_echs);

% Figure
% ------
close all
figure('Name','2D Example')
loglog(freqs_fft,sqrt(dsp_fft)*sqrt(ctn_noise_psd),'r');grid on;hold on;
loglog(freqs_ech,sqrt(dsp_ech),'b');
xl=get(gca,'XLim');set(gca,'XLim',[max(0.001,xl(1)) xl(2)])
yl=get(gca,'YLim');loglog(1/Te*[1 1],yl,'r--')
legend('PSD (1D-line cut from FFT2D)','PSD (on an axis of residuals)')
title({'';'\surd PSD (./\surdHz)'})
xlabel('Frequency (Hz)');ylabel('\surd PSD (./\surdHz)')