In general, we compute the power spectral density (PSD) of a noisy periodic signal. However, in this article, you will learn how to add noise directly to the PSD of a signal. This process is approximately equivalent to adding noise to a clean signal and then computing its PSD. Here, I will discuss both the theoretical background and the MATLAB implementation.
Steps
1. First, compute the Fast Fourier Transform (FFT) of the clean signal. Then, calculate the Power Spectral Density (PSD) from the FFT.
2. In our case, ensure that the PSD is in the linear scale. Next, compute the noise power from the given Signal-to-Noise Ratio (SNR) using:
noise_power = signal power / linear SNR
3. Then, generate Additive White Gaussian Noise (AWGN) using the formula:
AWGN noise = sqrt(noise_power) * randn
where randn generates a Gaussian-distributed signal with a mean of 0 and a variance of 1.
MATLAB Code
%% Define Parameters
fs = 1000; % Sampling frequency (Hz)
T = 0.2; % Time period of sine wave (s)
A = 1; % Amplitude
N = 1024; % Number of samples
t = linspace(-0.5, 0.5, N); % Time vector
f_sin = 5; % Frequency of sine wave (Hz)
%% Generate Periodic Sine Wave
sine_wave = A * sin(2 * pi * f_sin * t);
%% Compute PSD using FFT
Xf = fftshift(fft(sine_wave)); % Compute FFT and shift
PSD = abs(Xf).^2 / N; % Compute Power Spectral Density
%% Generate AWGN in Frequency Domain (Method 1)
snr_dB = 20; % SNR in dB
snr_linear = 10^(snr_dB/10); % Convert SNR to linear scale
signal_power = mean(PSD); % Approximate power of the original spectrum
noise_power = signal_power / snr_linear; % Compute noise power
noise_spectrum = sqrt(noise_power) .* (randn(size(PSD)) + 1j*randn(size(PSD))); % AWGN
%% Add AWGN Directly to PSD
noisy_PSD = PSD + abs(noise_spectrum).^2; % Add noise power to PSD
%% Generate AWGN in Time Domain (Method 2)
noise_time = sqrt(noise_power) * randn(size(sine_wave)); % AWGN in time domain
noisy_sine = sine_wave + noise_time; % Add noise to signal
%% Compute PSD of Noisy Sine Wave
Xf_noisy = fftshift(fft(noisy_sine)); % Compute FFT of noisy signal
PSD_noisy = abs(Xf_noisy).^2 / N; % Compute Power Spectral Density
%% Plot Results
freq = linspace(-fs/2, fs/2, N); % Frequency axis
figure;
% Plot Time-Domain Sine Wave
subplot(3,1,1);
plot(t, sine_wave, 'b', 'LineWidth', 1.5); hold on;
plot(t, noisy_sine, 'r', 'LineWidth', 1.2);
xlabel('Time (s)');
ylabel('Amplitude');
title('Sine Wave Before and After AWGN');
legend('Original Sine Wave', 'Noisy Sine Wave');
grid on;
% Plot PSD Comparison (Direct AWGN to PSD)
subplot(3,1,2);
plot(freq, 10*log10(PSD + eps), 'b', 'LineWidth', 1.5); hold on;
plot(freq, 10*log10(noisy_PSD + eps), 'r', 'LineWidth', 1.5);
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB)');
title('AWGN Added Directly to PSD');
legend('Original PSD', 'PSD with Direct AWGN');
grid on;
% Plot PSD Comparison (AWGN in Time Domain)
subplot(3,1,3);
plot(freq, 10*log10(PSD + eps), 'b', 'LineWidth', 1.5); hold on;
plot(freq, 10*log10(PSD_noisy + eps), 'g', 'LineWidth', 1.5);
xlabel('Frequency (Hz)');
ylabel('Power Spectral Density (dB)');
title('PSD: Original vs. PSD from Noisy Sine Wave');
legend('Original PSD', 'PSD from Noisy Signal');
grid on;
Output
Copy the MATLAB Code from here