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)}
Sampling frequency (fs): The rate at which the continuous-time signal is sampled (in Hz).
Total number of samples (N): The number of samples in the time-domain signal used for the DFT/FFT.
Suppose:
Sampling frequency = 1000 Hz
Number of samples = 500
Then:
Δf = 1000 / 500 = 2 Hz
This means the FFT result will contain frequency components spaced 2 Hz apart: 0 Hz, 2 Hz, 4 Hz, ..., up to fs.
- Increasing the number of samples (N) → improves frequency resolution
- Increasing the sampling frequency (fs) → worsens frequency resolution, but increases the total frequency range analyzed
MATLAB Script
% The code is written by SalimWireless.com
clear
close all
clc
fs = 1000; % 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*50 .* tt) + sin(2*pi*80 .* 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
% MINIMAL CHANGE 1: Multiply by standard deviation (sqrt of variance) for correct noise power
x = sqrt(noise_variance)*randn(L,1) + y; x = x(:); % sinusoid with additive Gaussian noise
% Plot results
figure
% Time-domain plot of the original signal
subplot(311)
% MINIMAL CHANGE 2: Swapped plot arguments and corrected labels
plot(tt, y,'r')
title('Original Message signal sin(2Ï€ * 50)t + sin(2Ï€ * 80)t (Time Domain)')
legend('Original signal')
xlabel('Time (s)')
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)')
web('https://www.salimwireless.com/search?q=psd%20fourier%20transform', '-browser');
% 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
Further Reading
- Bartlett Method for Spectral Estimation in MATLAB
- Periodogram for Spectral Estimation in MATLAB
- Welch Method for Spectral Estimation in MATLAB
- Calculation of SNR from FFT bins in MATLAB
- Add AWGN Directly to PSD in MATLAB
- How to Find the Fourier Transform of Any Signal
- Fourier Spectral Analysis
- Fourier Transform | Electronics Communication
- FFT Magnitude and Phase Spectrum using MATLAB
- Spectral Estimation Methods - Periodogram, Correlogram, Welch, Bartlett and Blackman-Tukey Methods