%
% file:      	noise_sounds.m, (c) Matthew Roughan, Sun Aug  1 2004
%

%%% first look at some white noise
fs = 44100;
T = 5;
x = (1:fs*T)/fs;
y = randn(size(x));

figure(10)
plot(x,y, 'linewidth', 3);
set(gca, 'xlim', [0 0.01]);
set(gca, 'linewidth', 3, 'fontsize', 18);
xlabel('time (seconds)');
print('-depsc', 'Plots/noise_white.eps');

z = fft(y);
w = abs(z).^2;
w = fftshift(w);
q = (-length(w)/2+1:length(w)/2)/T;

figure(1)
semilogy(q/1000, w, 'linewidth', 3);
set(gca, 'linewidth', 3);
set(gca, 'fontsize', 18);
set(gca, 'xlim', [0 fs/2000]);
xlabel('frequency (kHz)');
print('-depsc', 'Plots/noise_white_fft.eps');

wavwrite(y, fs, 'Plots/noise_white.wav');







%%% now filter the noise with a simple MA
ma = ones(11,1)/sqrt(11);
y1 = filter(ma, 1, y);
f1 = zeros(size(y1));
f1(1:length(ma)) = ma;
figure(20)
plot(x,y1, 'linewidth', 3);
set(gca, 'xlim', [0 0.01]);
set(gca, 'linewidth', 3);
set(gca, 'fontsize', 18);
xlabel('time (seconds)');
print('-depsc', 'Plots/noise_ma.eps');

z1 = fft(y1);
z2 = fft(f1);
w1 = abs(z1).^2;
w1 = fftshift(w1);
w2 = abs(z2).^2;
w2 = fftshift(w2);
w2 = w2 * (max(w1)/max(w2))/10;
q = (-length(w)/2+1:length(w)/2)/T;

figure(2)
hold off
semilogy(q/1000, w1, 'linewidth', 3);
hold on
semilogy(q/1000, w2, 'r', 'linewidth', 2);
set(gca, 'linewidth', 3);
set(gca, 'fontsize', 18);
set(gca, 'xlim', [0 fs/2000]);
xlabel('frequency (kHz)');
print('-depsc', 'Plots/noise_ma_fft.eps');

wavwrite(y1, fs, 'Plots/noise_ma.wav');


%%% now filter the noise with a longer MA
ma = ones(111,1)/sqrt(111);
y1 = filter(ma, 1, y);
f1 = zeros(size(y1));
f1(1:length(ma)) = ma;
figure(30)
plot(x,y1, 'linewidth', 3);
set(gca, 'xlim', [0 0.01]);
set(gca, 'linewidth', 3);
set(gca, 'fontsize', 18);
xlabel('time (seconds)');
print('-depsc', 'Plots/noise_ma_2.eps');

z1 = fft(y1);
z2 = fft(f1);
w1 = abs(z1).^2;
w1 = fftshift(w1);
w2 = abs(z2).^2;
w2 = fftshift(w2);
w2 = w2 * (max(w1)/max(w2))/10;
q = (-length(w)/2+1:length(w)/2)/T;

figure(3)
hold off
semilogy(q/1000, w1, 'linewidth', 3);
hold on
semilogy(q/1000, w2, 'r', 'linewidth', 2);
set(gca, 'linewidth', 3);
set(gca, 'fontsize', 18);
set(gca, 'xlim', [0 fs/2000], 'ylim', 10.^[0 10]);
xlabel('frequency (kHz)');
print('-depsc', 'Plots/noise_ma_fft_2.eps');

wavwrite(y1, fs, 'Plots/noise_ma_2.wav');

%%% now filter a difference
difference = [1 -1];
y1 = filter(ma, 1, y);
f1 = zeros(size(y1));
f1(1:length(difference)) = difference;
y1 = filter(difference, 1, y);
figure(40)
plot(x,y1, 'linewidth', 3);
set(gca, 'xlim', [0 0.01]);
set(gca, 'linewidth', 3);
set(gca, 'fontsize', 18);
xlabel('time (seconds)');
print('-depsc', 'Plots/noise_diff.eps');


z1 = fft(y1);
z2 = fft(f1);
w1 = abs(z1).^2;
w1 = fftshift(w1);
w2 = abs(z2).^2;
w2 = fftshift(w2);
w2 = w2 * (max(w1)/max(w2))/10;
q = (-length(w)/2+1:length(w)/2)/T;

figure(4)
hold off
semilogy(q/1000, w1, 'linewidth', 3);
hold on
semilogy(q/1000, w2, 'r', 'linewidth', 2);
set(gca, 'linewidth', 3);
set(gca, 'fontsize', 18);
set(gca, 'xlim', [0 fs/2000], 'ylim', 10.^[0 10]);
xlabel('frequency (kHz)');
print('-depsc', 'Plots/noise_diff_fft.eps');

wavwrite(y1, fs, 'Plots/noise_diff.wav');



