% Pink Noise Maker
% Steven Backer
% MUS 220a HW #5
% November 2004

% Filter Coefficients Adapted from
%	Smith, Julius O. Introduction to Digital Filters, May 2004 Draft,
%	Center for Computer Research in Music and Acoustics (CCRMA), 
%	Stanford University.

clear all;

% Filter Coefficients
q = [0.98223142883799 0.83265674642178 0.10798087989066];
p = [0.99516895763438 0.94384178587182 0.55594525865504];
B = zeros(1,4);
A = zeros(1,4);
B = poly(q);
A = poly(p);

NoiseL = 120;  % number of samples to synthesize

sopranoNoise = randn(1,NoiseL);
altoNoise = randn(1,NoiseL);
tenorNoise = randn(1,NoiseL);

pinkSoprano = filter(B,A,sopranoNoise);
pinkAlto = filter(B,A,altoNoise);
pinkTenor = filter(B,A,tenorNoise);

save psn360 pinkSoprano -ascii 
save pan360 pinkAlto -ascii 
save ptn360 pinkTenor -ascii 

% analyze_1f
fs = 48000; % If this were "noise"
n=[1:fs/2];
ideal = (1./(n));
refdB = 0.0004;  % arbitrary
semilogx(n, 20*log10(ideal/refdB),'m')
grid off
hold on

[Pss, w] = pwelch(pinkSoprano);
semilogx(w*fs/(2*pi), 20*log10(Pss), 'r')
Paa = pwelch(pinkAlto);
semilogx(w*fs/(2*pi), 20*log10(Paa), 'g')
Ptt = pwelch(pinkTenor);
semilogx(w*fs/(2*pi), 20*log10(Ptt))

legend('1/f','Soprano','Alto','Tenor')
axis([w(2)*fs/(2*pi) 24000 -30 30])
hold off
xlabel('Frequency [Hz]')
ylabel('PSD [dB]')
title('Estimated Power Spectal Density of Generated Pink Noise')

% save session
save
