%% Two-Channel Brownian Noise
% 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 two channels of Brownian noise by setting Color to 'brown' and NumChannels to 2.

cn = dsp.ColoredNoise('brown','SamplesPerFrame',2048,...
    'NumChannels',2);
x = cn();
subplot(2,1,1)
plot(x(:,1)); title('Channel 1'); axis tight;
subplot(2,1,2)
plot(x(:,2)); title('Channel 2'); axis tight;

%
% The sampling frequency is 1 Hz. Obtain Welch PSD estimates for both channels. The fourth argument of pwelch, NFFT, which is the number of FFT points, is empty. Hence, NFFT is set to 256. For even NFFT, The number of FFT points used to calculate the PSD estimate is (NFFT/2+1), which equals 129.

Fs = 1;
Pxx = zeros(129,size(x,2));
for nn = 1:size(x,2)
[Pxx(:,nn),F] = pwelch(x(:,nn),hamming(128),[],[],Fs,'psd');
end

%
% Construct the theoretical PSD of a Brownian process. Plot the theoretical PSD along with both realizations on a log-log plot. Use a base-2 logarithmic scale for the frequency axis and plot the power spectral densities in dB.

PSDBrownian = 1./F(2:end).^2;
figure;
plot(log2(F(2:end)),10*log10(PSDBrownian),'k-.','linewidth',2);
hold on;
plot(log2(F(2:end)),10*log10(Pxx(2:end,:)));
xlabel('log_2(Hz)'); ylabel('dB');
grid on;
legend('Theoretical PSD','Channel 1', 'Channel 2');

% Copyright 2012 The MathWorks, Inc.