%% PSD of Pink Noise Realization
% Note: If you are using R2016a or an earlier release, replace each call to the object with the equivalent step syntax. For example, obj(x) becomes step(obj,x).
% Generate a pink noise signal 2048 samples in length. The sampling frequency is 1 Hz. Obtain an estimate of the power spectral density using Welch's overlapped segment averaging.

cn = dsp.ColoredNoise('pink','SamplesPerFrame',2048);
x = cn();
Fs = 1;
[Pxx,F] = pwelch(x,hamming(128),[],[],Fs,'psd');

%
% Construct the theoretical PSD of the pink noise process.

PSDPink = 1./F(2:end);

%
% Display the Welch PSD estimate of the noise along with the theoretical PSD on a log-log plot. Plot the frequency axis  with a base-2 logarithmic scale to clearly show the octaves. Plot the PSD estimate in dB, 10\log_{10}.

plot(log2(F(2:end)),10*log10(Pxx(2:end)))
hold on
plot(log2(F(2:end)),10*log10(PSDPink),'r','linewidth',2)
xlabel('log_2(Hz)')
ylabel('dB')
title('Pink Noise')
grid on
legend('PSD estimate','Theoretical pink noise PSD')
hold off

% Copyright 2012 The MathWorks, Inc.