function NoiseColorTest
% Noise Color Test, version 3 (January 2017), by toh@umd.edu
% The frequency distribution of noise, designated by noise color,
% substantially effects the ability of smoothing to reduce noise. The
% Matlab/Octave function “NoiseColorTest.m” compares the effect of a
% 20-point boxcar (unweighted sliding average) smooth on the standard
% deviation of white, pink, and blue noise, all of which have an original
% unsmoothed standard deviation of 1.0. Because smoothing is a low-pass
% filter process, it effects low frequency (pink and red) noise less, and
% effects high-frequency (blue and violet) noise more, than it does white
% noise. Figure 1 shows the five type of noises before smoothing, figure 3
% shown their frequency spectra, and figure 2 shows the five type of noises
% after smoothing. 
format compact
% User-variable constants:
x=1:1000;
Noise=1;
SmoothWidth=20;
SmoothType=1; % 1=rectangular (sliding average); 2=triangle; 3=Gaussian

WhiteNoiseArray=Noise*(randn(size(x))); % normal white constant noise
BlueNoiseArray=Noise*bluenoise(length(x)); % Blue noise
PinkNoiseArray=Noise*pinknoise(length(x)); % Pink (1/f) noise
RedNoiseArray=cumsum(WhiteNoiseArray)./(std(cumsum(WhiteNoiseArray)));
VioletNoiseArray=Noise*violetnoise(length(x)-1); % Violet noise
size(VioletNoiseArray)

SmoothedWhiteNoise=fastsmooth(WhiteNoiseArray,SmoothWidth,SmoothType);
SmoothedPinkNoise=fastsmooth(PinkNoiseArray,SmoothWidth,SmoothType);
SmoothedBlueNoise=fastsmooth(BlueNoiseArray,SmoothWidth,SmoothType);
SmoothedRedNoise=fastsmooth(RedNoiseArray,SmoothWidth,SmoothType);
SmoothedVioletNoise=fastsmooth(VioletNoiseArray,SmoothWidth,SmoothType);

disp('Standard Deviations:')
BeforeSmoothingSTD=stdev(WhiteNoiseArray)
SmoothedWhiteSTD=stdev(SmoothedWhiteNoise)
SmoothedPinkSTD=stdev(SmoothedPinkNoise)
SmoothedBlueSTD=stdev(SmoothedBlueNoise)
SmoothedRedSTD=stdev(SmoothedRedNoise)
SmoothedVioletSTD=stdev(SmoothedVioletNoise)

figure(1)
clf
subplot(2,3,2)
plot(x,WhiteNoiseArray)
axis([0 length(x) -3 3])
title('White Noise')
xlabel(['Standard deviation = ' num2str(stdev(WhiteNoiseArray))] )
subplot(2,3,3)
plot(x,PinkNoiseArray)
axis([0 length(x) -3 3])
title('Pink (1/f) Noise')
xlabel(['Standard deviation = ' num2str(stdev(WhiteNoiseArray))] )
subplot(2,3,4)
plot(x,BlueNoiseArray)
axis([0 length(x) -3 3])
title('Blue Noise')
xlabel(['Standard deviation = ' num2str(stdev(PinkNoiseArray))] )
subplot(2,3,5)
plot(x,RedNoiseArray)
title('Red Noise')
xlabel(['Standard deviation = ' num2str(stdev(RedNoiseArray))] )
axis([0 length(x) -3 3])
subplot(2,3,6)
plot(x,VioletNoiseArray)
title('Violet Noise')
xlabel(['Standard deviation = ' num2str(stdev(VioletNoiseArray))] )
axis([0 length(x) -3 3])

figure(2)
subplot(2,3,1)
plot(x,WhiteNoiseArray,'c')
hold on
switch SmoothType,
    case 1
        plot(x,rectanglepulse(x,length(x)/2,SmoothWidth),'r')
    case 2
        plot(x,triangle(x,length(x)/2,SmoothWidth),'r')
    case 3
        plot(x,gaussian(x,length(x)/2,SmoothWidth),'r')
end
axis([0 length(x) -3 3])
hold off
title('White noise (Red=smoothing funciton)')
xlabel(['Width = '  num2str(SmoothWidth)    '   Type = '  num2str(SmoothType) ] )
subplot(2,3,2)
plot(x,SmoothedWhiteNoise,'c')
axis([0 length(x) -3 3])
title('Smoothed White Noise')
xlabel(['Standard deviation = ' num2str(SmoothedWhiteSTD)] )
subplot(2,3,3)
plot(x,SmoothedPinkNoise,'c')
axis([0 length(x) -3 3])
title('Smoothed Pink (1/f) Noise')
xlabel(['Standard deviation = ' num2str(SmoothedPinkSTD)] )
subplot(2,3,4)
plot(x,SmoothedBlueNoise,'c')
axis([0 length(x) -3 3])
title('Smoothed Blue Noise')
xlabel(['Standard deviation = ' num2str(SmoothedBlueSTD)] )
subplot(2,3,5)
plot(x,SmoothedRedNoise,'c')
title('Smoothed Red Noise')
xlabel(['Standard deviation = ' num2str(SmoothedRedSTD)] )
axis([0 length(x) -3 3])
subplot(2,3,6)
plot(x,SmoothedVioletNoise,'c')
title('Smoothed Violet Noise')
xlabel(['Standard deviation = ' num2str(SmoothedVioletSTD)] )
axis([0 length(x) -3 3])

figure(3);
clf
subplot(2,3,2)
PlotFrequencySpectrum(x,WhiteNoiseArray',4,0,0);
title('WhiteNoiseArray')
xlabel('Log-log frequency spectrum')
subplot(2,3,3)
PlotFrequencySpectrum(x,PinkNoiseArray',4,0,0);
title('PinkNoiseArray')
xlabel('Log-log frequency spectrum')
subplot(2,3,4)
PlotFrequencySpectrum(x,BlueNoiseArray',4,0,0);
title('BlueNoiseArray')
xlabel('Log-log frequency spectrum')
subplot(2,3,5)
PlotFrequencySpectrum(x,RedNoiseArray',4,0,0);
title('RedNoiseArray')
xlabel('Log-log frequency spectrum')
subplot(2,3,6)
PlotFrequencySpectrum(x,VioletNoiseArray',4,0,0);
title('VioletNoiseArray')
xlabel('Log-log frequency spectrum')

function y=whitenoise(x)
% Random noise with white power spectrum with mean zero 
% and unit standard deviation, equal in length to x
% Tom O'Haver, 2008
% Example:
% model=gaussian([1:10deriv00],500,200);
% model=model+.1.*whitenoise(model);
% plot(model)
% 
y=randn(size(x));

function by=bluenoise(n)
% Random noise with blue power spectrum with mean zero 
% and unit standard deviation. n is number of points.
% Tom O'Haver, 2008
% Example:
% model=gaussian([1:1000],500,200);
% model=model+.1.*bluenoise(length(model));
% plot(model)
x=1:n;
y=randn(size(x));  % Random normally-distributed white noise
ry=deriv1(y);
by=(ry./stdev(ry)); % Normalize to unit standard deviation

function stddev=stdev(a)
% Octave and Matlab compatible standard deviation function
sa=size(a);
reshape(a,1,length(a));
if sa(1)>sa(2),
  stddev=std(a);
else
  stddev=(std(a'));
end;

function ry=pinknoise(n)
% Random noise with pink (1/f) power spectrum with mean zero 
% and unit standard deviation. n is number of points.
% Tom O'Haver, 2008
% Example:
% model=gaussian([1:1000],500,200);
% model=model+.1.*pinknoise(length(model));
% plot(model)
x=[1:n];
y=randn(size(x));  % Random normally-distributed white noise
% Fourier filter 
fy=fft(y); % Compute Fourier transform of signal y
% Compute filter shape
lft1=[1:(length(fy)/2)+1];
lft2=[(length(fy)/2):length(fy)];
ffilter1=ones(size(lft1))./(sqrt(lft1));
ffilter2=ones(size(lft2))./(sqrt(lft2));
ffilter=[ffilter1,ffilter2];
if length(fy)>length(ffilter), ffilter=[ffilter ffilter(1)];end
ffy=fy.*ffilter(1:length(fy));  % Multiply filter by Fourier transform of signal
ry=real(ifft(ffy)); % Inverse transform to recover filtered signal 'ry'
ry=((ry-mean(ry))./stdev(ry)); % Normalize to unit standard deviation

function vy = violetnoise(x)
% Random noise with high-frequency weighted power spectrum with mean zero
% and unit standard deviation. Length equal to the length of x
y=0:x;
for n=1:2:length(y)-1,
    rn=abs(randn());
    y(n)=rn;
    y(n+1)=-rn;
end
vy=(y-mean(y))./stdev(y); % Normalize to unit standard deviation

function d=deriv1(a)
% First derivative of vector using simple adjacent differences.
% Example: deriv([1 1 1 2 3 4]) yeilds [0 0 0 1 1 1]
%  T. C. O'Haver, 2013.
n=length(a);
d=zeros(size(a));
d(1)=a(2)-a(1);
for j = 2:n;
  d(j)=a(j)-a(j-1);
end
