Discrete Power Spectral Density

From Wikimization

%Demonstrate receding noise floor in digital domain for fixed analog power spectral density.
%For COA paper.  powerSpectralDensity.m
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 * 22400);  % RMS thermal noise voltage from one-sided spectrum

% 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/1000), 'k']);
end
% Labels & legend
xlabel('Frequency [Hz]');
ylabel('1k\Omega Discrete Power Spectral Density [dBFS]');
title ('Thermal Noise Spectral Density for Decade Record Lengths');
% Create legend and retrieve its icon handles
[lgd,icons] = legend('show','Location','southwest');
% icons is an array―lines for curves, patches for fills, etc.
% Thicken only the line icons; leave markers and patches alone
for ii = 1:numel(icons)
    if strcmp(icons(ii).Type,'line')          % icon is a line sample
        icons(ii).LineWidth = 2.5;            % make it thicker
    end
end
lgd.Units    = 'normalized';                   % figure-relative coordinates
lgd.Position = [0.4275 0.18 0.18 0.12];        % [left  bottom  width  height]
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'])