Discrete Power Spectral Density

From Wikimization

%Demonstrate receding noise floor in digital domain for fixed analog power spectral density
clearvars; clc; close all
disp('This script demonstrates how thermal noise spectral level') 
disp('appears to fall in the frequency domain as record length increases.');
Fs = 48000;  % Sample rate
R = 1000;    % Resistance [Ohms]
k = 1.380649e-23;  % Boltzmann constant [J/K]
T_R = 300;   % Resistor Temperature 300°K roughly equivalent to 27°C
Vn_rms = sqrt(4 * k * T_R * R * Fs/2);  % RMS thermal noise voltage

% Record lengths
N_values = [Fs, 10*Fs, 100*Fs, 1000*Fs];

figure; hold on;
% Loop over different record lengths
for i = 1:numel(N_values)
   N = N_values(i);       %Number of samples for current record.
   freq = Fs*(0:N-1)'/N;  %frequency vector

   % Generate thermal noise signal (standard Gaussian white noise)
   noise = Vn_rms * randn(N, 1);
   Y = fft(noise)/N;   %no window

   %Compute average spectral level
   DFT_square            = real(Y.*conj(Y)); 
   Avg_Spectral_Level(i) = 10*log10(sum(DFT_square)/N);  %Average Spectral Level recedes by 10dB at each pass.
   Noise_Level(i)        = 10*log10(sum(DFT_square));    %Absolute Noise Level is fixed and independent of N.
   plot(freq,              10*log10(    DFT_square), 'DisplayName', ...
                                                    ['Avg Spectral Level: ',      num2str(Avg_Spectral_Level(i),'%4.0f'), 'dB' ... 
                                                     ',  Absolute Noise Level: ', num2str(Noise_Level(i),       '%4.0f'), 'dB' ...
                                                     ',  N: ',                    num2str(N)]);
end
% Labels & legend
xlabel('Frequency [Hz]');
ylabel('Discrete Power Spectral Density [dB/Hz]');
title ('Thermal Noise Spectral Density for Decade Record Lengths');
legend('show', 'Location', 'southwest');
xticks(0:9600:Fs);
ax = gca;
ax.XAxis.Exponent = 0;  %disable scientific notation along x axis
ax.XTickLabel = arrayfun(@num2str, ax.XTick, 'UniformOutput', false);
xlim([0, Fs]);
ylim([-210,-160]);
grid on; hold off;

disp(' ')
disp('Sample rate is fixed at 48 kHz.')  
disp('Record lengths are Fs, 10Fs, 100Fs, and 1000Fs.');
disp(' ')
disp(['Noise level w.r.t unit amplitude sinusoid: ' num2str(20*log10(10^(Noise_Level(i)/20)/(1/sqrt(2))),'%4.0f') 'dB'])