Skip to main content

Power Spectral Density Calculation Using FFT in MATLAB


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

  1. Firstly, calculate the fast Fourier transform (FFT) of a signal.
  2. Then, calculate the Fourier magnitude (absolute value) of the signal.
  3. Square the Fourier magnitude to get the power spectrum.
  4. To calculate the Power Spectral Density (PSD), divide the squared magnitude by the product of the sampling frequency (fs) and the total number of samples (N).
    Formula: PSD = |FFT|^2 / (fs * N)

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 the frequency resolution is:
Δ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 (Nyquist limit)

MATLAB Script

matlab_script.m

% 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

% 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)
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 Periodogram 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

Power Spectral Density Periodogram output plot

Power Spectral Density output visualization

Interactive PSD Simulator

Experiment with signal parameters and visualize the Power Spectral Density in real-time below.

Parameters


Base Signal



Input Signal

Magnitude Plot

Phase Plot



10
rectangular


Contact Us

Name

Email *

Message *

Popular Posts

Q-function in BER vs SNR Calculation (with Simulation)

Q-function in BER vs. SNR Calculation In digital communications and signal processing, the Q-function plays a significant role in predicting system reliability. It allows engineers to quantify the probability that Gaussian noise will exceed a specific threshold, causing a bit error. What is the Q-function? The Q-function is a mathematical function representing the tail probability of the standard normal (Gaussian) distribution. It is the complementary cumulative distribution function (CCDF) of a standard Gaussian distribution. Q(x) = (1 / √(2Ï€)) ∫â‚“∞ e^(-t² / 2) dt Q-Function Interactive Simulator Move the slider to see how the "Tail Probability" (the area in red) changes. This area represents the Probability of Error (BER) . Threshold Distance ( x ) — (Simulates Increasing SNR) x = 1.0 Q(x) = 0.1587 ...

BER vs SNR for M-ary QAM, M-ary PSK, QPSK, BPSK, ...(MATLAB Code + Simulator)

Bit Error Rate (BER) & SNR Guide Analyze communication system performance with our interactive simulators and MATLAB tools. 📘 Theory 🧮 Simulators 💻 MATLAB Code 📚 Resources BER Definition SNR Formula BER Calculator MATLAB Comparison 📂 Explore M-ary QAM, PSK, and QPSK Topics ▼ 🧮 Constellation Simulator: M-ary QAM 🧮 Constellation Simulator: M-ary PSK 🧮 BER calculation for ASK, FSK, and PSK 🧮 Approaches to BER vs SNR What is Bit Error Rate (BER)? The BER indicates how many corrupted bits are received compared to the total number of bits sent. It is the primary figure of merit f...

Frequency Shift Keying (FSK) Modulation & Demodulation (with Simulation)

Frequency Shift Keying (FSK) Theoretical Foundations: Frequency Shift Keying (FSK) is a discrete frequency modulation scheme wherein the digital information is encoded via instantaneous shifts in the carrier signal's frequency. The fundamental implementation is Binary FSK (BFSK), which maps binary data onto two distinct, discrete spectral states. A binary '1' (the "mark" state) is represented by a carrier frequency \( f_1 \), while a binary '0' (the "space" state) corresponds to frequency \( f_2 \). Each symbol is sustained for a bit interval denoted by \( T_b \). FSK Transmitter Characterization: The mathematical model for the modulated BFSK output \( s(t) \) is defined as: \[ s(t) = \begin{cases} A_c \cos(2\pi f_1 t), & \text{for } m = 1 \\ A_c \cos(2\pi f_2 t), & \text{for } m = 0 \end{cases} \] ...

RMS Delay Spread, Excess Delay Spread and Multi-path ...(with MATLAB + Simulator)

📘 Overview of Delay Spread and Multi-path 🧮 Excess Delay spread 🧮 Power delay Profile 🧮 RMS Delay Spread 📚 Further Reading 📂 Other Topics on RMS Delay Spread, Excess Delay ... 🧮 Multipath Components or MPCs 🧮 Online Simulator for Calculating RMS Delay Spread 🧮 Why is there significant multipath in the case of very high frequencies? 🧮 Why RMS Delay Spread is essential for wireless communication? 🧮 Why the Power Delay Profile is essential? 🧮 MATLAB Codes for Calculating Different Types of delay Spreads Delay Spread, Excess Delay Spread, and Multipath (MPCs) The fundamental distinction between wireless and wired connections is that in wireless connections signal reaches at receiver thru multipath signal propagation rather than directed transmission like co-axial cable. Wireless Communication has no set communication path between the transmitter and the receiver. The line...

FFT Butterfly Method Explained (with Example of 4-point DFT)

  FFT Using Butterfly Method Given: x[n] = {0, 1, 2, 3} Step 1: Split into Even & Odd Even indices: x e = {0, 2} Odd indices: x o = {1, 3} Step 2: 2-point DFT For any {a, b}: DFT = {a + b, a - b} Even Part: E = {0+2, 0-2} = {2, -2} Odd Part: O = {1+3, 1-3} = {4, -2} Step 3: Combine Using Butterfly X[k] = E[k] + W k O[k] X[k + N/2] = E[k] - W k O[k] For N = 4: W 0 = 1 W 1 = -j Final Calculations X[0] = 2 + 4 = 6 X[2] = 2 - 4 = -2 X[1] = -2 + (-j)(-2) = -2 + 2j X[3] = -2 - (-j)(-2) = -2 - 2j Final Answer: X[k] = {6, -2 + 2j, -2, -2 - 2j} Try Interactive Online Simulations Interactive FFT Online Simulator (For understanding Fundamentals)  Interactive FFT Online Simulator (Analyze .CSV, .MP3, .MP4, etc. Further Reading Fourier Transform OFDM Return to Fourier Transform Main Page →

AM Modulation Online Simulator

Amplitude Modulation Simulator s AM (t) = A c [1 + k a m(t)] cos(ω c t) where, ω = 2πf & k a = Amplitude Sensitivity Modulation index, μ = k a A m Message Frequency (fm): Carrier Frequency (fc): Carrier Amplitude (Ac): Modulation Index (m = Am / Ac):

Pulse Width Modulation (PWM)

Pulse-width modulation (PWM), or pulse-duration modulation (PDM), is a method of controlling the average power delivered by an electrical signal.   Fig: An example of PWM in an idealized inductor driven by a blue line voltage source modulated as a series of sawtooth pulses, resulting in a red line current in the inductor.    Generating a PWM Signal The simplest way to generate a PWM signal is the intersection method, which requires only a sawtooth or a triangle waveform (easily generated using a simple oscillator) and a comparator. When the value of the reference signal is more than the modulation waveform, the PWM signal (magenta) is in the high state; otherwise, it is in the low state.      Duty cycle A low duty cycle equates to low power because the power is off for most of the time; the word duty cycle reflects the ratio of "on" time to the regular interval or "period" of time. The duty cycle is measured in percent, with 100% representing full o...