Power spectral density (PSD) tells us how the power of a signal is distributed across different frequency components, whereas Fourier Magnitude gives you the amplitude (or strength) of each frequency component in the signal.
Steps to calculate the PSD of a signal
- Firstly, calculate the first Fourier transform (FFT) of a signal
- Then, calculate the Fourier magnitude of the signal
- The power spectrum is the square of the Fourier magnitude
- To calculate power spectrum density (PSD), divide the power spectrum by the total number of samples and the frequency resolution. {Frequency resolution = (sampling frequency / total number of samples)}
MATLAB Script
% The code is written by SalimWireless.com
clearclose all
clc
fs = 40000; % sampling frequency
T = 1; % total recording time
L = T .* fs; % signal length
tt = (0:L-1)/fs; % time vector
ff = (0:L-1)*fs/L;
y = sin(2*pi*100 .* tt) .* sin(2*pi*1000 .* tt); y = y(:); % reference sinusoid
% Allow user to input SNR in dB
snr_db = input('Enter the SNR (in dB): '); % User input for SNR
snr_linear = 10^(snr_db / 10); % Convert SNR from dB to linear scale
% Calculate noise variance based on SNR
signal_power = mean(y.^2); % Calculate signal power
noise_variance = signal_power / snr_linear; % Calculate noise variance
x = noise_variance*randn(L,1) + y; x = x(:); % sinusoid with additive Gaussian noise
% Plot results
figure
% Manual Power Spectral Density plots
subplot(311)
plot(y, tt,'r')
title('Original Message signal')
legend('Original signal')
xlabel('Frequency (Hz)')
ylabel('Amplitude')
% Manual Power Spectral Density plots
subplot(312)
[psd_y, f_y] = manualPSD(y, fs); % PSD of the original signal
plot(f_y,10*log10(psd_y),'r')
title('Power Spectral Density')
legend('Original signal PSD')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
% Manual Power Spectral Density plots
subplot(313)
[psd_x, f_x] = manualPSD(x, fs); % PSD of the noisy signal
plot(f_x,10*log10(psd_x),'k')
title('Power Spectral Density')
legend('Noisy signal PSD')
xlabel('Frequency (Hz)')
ylabel('Power/Frequency (dB/Hz)')
% Manual PSD calculation function
function [psd, f] = manualPSD(signal, fs)
N = length(signal); % Signal length
fft_signal = fft(signal); % FFT of the signal
fft_signal = fft_signal(1:N/2+1); % Take only the positive frequencies
psd = (1/(fs*N)) * abs(fft_signal).^2; % Compute the power spectral density
psd(2:end-1) = 2*psd(2:end-1); % Adjust the PSD for the one-sided spectrum
f = (0:(N/2))*fs/N; % Frequency vector
end
Output
