%% Octave Smoothing of White and Pink Noise
% Generate {10}^5 samples of white Gaussian noise sampled at 44.1 kHz. Create a signal of pink noise by filtering the white noise with a filter whose zeros and poles are all on the positive x-axis.

N = 1e5;
fs = 44.1e3;
wn = randn(N,1);

z = [0.982231570015379 0.832656605953720 0.107980893771348]';
p = [0.995168968915815 0.943841773712820 0.555945259371364]';
[b,a] = zp2tf(z,p,1);

pn = filter(b,a,wn);

%
% Compute the Welch estimate of the power spectral density for both signals. Divide the signals into 2048-sample segments, specify 50% overlap between adjoining segments, window each segment with a Hamming window, and use 4096 DFT points.

[pxx,f] = pwelch([wn pn],hamming(2048),1024,4096,fs);

% Display the spectral densities over a frequency band ranging from 200 Hz to the Nyquist frequency. Use a logarithmic scale for the frequency axis.

pwelch([wn pn],hamming(2048),1024,4096,fs)
ax = gca;
ax.XScale = 'log';
xlim([200 fs/2]/1000)
legend('White','Pink')

%
% Compute and display the octave spectra of the signals. Use the same frequency range as in the previous plot. Specify six bands per octave and compute the spectra using 8th-order filters.

poctave(pxx,fs,f,'BandsPerOctave',6,'FilterOrder',8,'FrequencyLimits',[200 fs/2],'psd')
legend('White','Pink')

% Copyright 2012 The MathWorks, Inc.