
% Adam Wang, 17 Mar 2008

% Makes colored noise and their NPS


randn('state',2);

imSize = 256;

white = randn(imSize,imSize);
White = fft2c(white);

x = ones(imSize,1)*(1:imSize);
y = x';
dist = sqrt((x-(imSize/2+1)).^2 + (y-(imSize/2+1)).^2); % distance from center of image

sqrtRamp = sqrt(dist/max(dist(:))); % normalized 1/sqrt(rho) filter
redFilt = sqrt(1-dist/max(dist(:)));

B = White.*sqrtRamp;
b = ifft2c(B);

blue = b*sqrt(mean(white(:).^2)/mean(b(:).^2)); % scale blue noise to same MSE as white
Blue = fft2c(blue);

R = White.*redFilt;
r = ifft2c(R);

red = r*sqrt(mean(white(:).^2)/mean(r(:).^2)); % scale red noise to same MSE as white
Red = fft2c(red);

Max = max(max([white red blue]));
Min = min(min([white red blue]));

imwrite(uint8((white-Min)/(Max-Min)*256), 'C:\Documents and Settings\Adam Wang\My Documents\_Second Year\PSYCH221\Project\data\white.png', 'png');
imwrite(uint8((red-Min)/(Max-Min)*256), 'C:\Documents and Settings\Adam Wang\My Documents\_Second Year\PSYCH221\Project\data\red.png', 'png');
imwrite(uint8((blue-Min)/(Max-Min)*256), 'C:\Documents and Settings\Adam Wang\My Documents\_Second Year\PSYCH221\Project\data\blue.png', 'png');


MSE_white = mean(white(:).^2);
MSE_blue = mean(blue(:).^2);
MSE_red = mean(red(:).^2);

NPSwhite = zeros(1,imSize/2);
NPSblue = zeros(1,imSize/2);

for j = 1:16
    w = randn(imSize,imSize);
    W = fft2c(w);
    B = W.*sqrtRamp;
    b = ifft2c(B);
    bb = b*sqrt(mean(w(:).^2)/mean(b(:).^2));
    BB = fft2c(bb);
    
    for i = 1:imSize/2
        NPSwhite(i) = NPSwhite(i) + mean(abs(W(and(dist>=i-1,dist<i))).^2);
        NPSblue(i) = NPSblue(i) + mean(abs(BB(and(dist>=i-1,dist<i))).^2);
    end
end

NPSwhite = NPSwhite/imSize^2/16; %avg NPS of 16 realizations, function of rho
NPSblue = NPSblue/imSize^2/16;  %rho is radial frequency

