Poisson noise matlab codes

% Program to improve S/N ratio in small factor steps of 2^0.1 = 1.0718 of a

% 16-bit tif photograph (copyright Nyah Willmott)

clear; close all;

vid = VideoWriter('poissonNoise.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 15;

open(vid);

set(0,'defaultfigurecolor',[1 1 1]);

set(gca,'linewidth',7);

I = imread('zurich.tif'); % input image

J = rgb2gray(I); % convert to b&w

K = imresize(J,0.35); % rescale image to more practical size

for i = 16:-0.1:0

factor = 2^i;

L = K/factor; % reduce bit range. Reduces each pixel to round(val/factor)

M = imnoise(L,'poisson'); % Add Poisson noise

N = M*factor; % Multiply up again, now with noise added

imshow(N);

bits = round(16.4 - i);

str1 = num2str(bits,'% 2i');

if (bits >=1 )

if (bits == 1)

str2 = ' bit';

else

str2 = ' bits';

end

strTot = [str1,str2];

hText = text(1600, 1190, strTot, 'FontSize',16,'Color','white');

end

frame = getframe(gcf);

writeVideo(vid,frame);

end

% Output the movie as an mpg file

close(vid);

% Program to improve S/N ratio in small factor steps of 2^0.1 = 1.0718 of a

% 16-bit tif photograph (copyright Nyah Willmott)

clear; close all;

vid = VideoWriter('poissonNoiseCroppedRegion.mp4','MPEG-4');

vid.Quality = 100;

vid.FrameRate = 15;

open(vid);

set(0,'defaultfigurecolor',[1 1 1]);

set(gca,'linewidth',7);

I = imread('zurich.tif','PixelRegion',{[1554,1680],[1760,1873]}); % input image

J = rgb2gray(I); % convert to b&w

for i = 16:-0.1:0

factor = 2^i;

L = J/factor; % reduce bit range. Reduces each pixel to round(val/factor)

M = imnoise(L,'poisson');

N = M*factor;

imshow(N);

bits = round(16.4 - i);

str1 = num2str(bits,'% 2i');

if (bits >=1 )

if (bits == 1)

str2 = ' bit';

else

str2 = ' bits';

end

strTot = [str1,str2];

hText = text(1600, 1190, strTot, 'FontSize',16,'Color','white');

end

frame = getframe(gcf);

writeVideo(vid,frame);

end

% Output the movie as an mpg file

close(vid);

% Program to generate signal with noise and background with noise, fit them

% via a least squares fitting function provided by matlab, and plot the

% progress of the fit with number of scans.

% The function is a constant background with amplitude = 25, a peak signal

% with height = 2, centred at x = 250 and with a SD value of 16.

clear;

close all;

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');

vid = VideoWriter('snr.mp4','MPEG-4');

vid.FrameRate = 60; % Default 30

vid.Quality = 100; % Default 75

open(vid);

set(0,'defaultfigurecolor',[1 1 1]);

set(gca,'linewidth',7);

philBlue = [0.3 0.28 0.55];

numScans = 2000; % Number of scans

scanLength = 500; % Length of any given scan

sigPos = scanLength/2; % Place peak signal ideally at centre of scan

sigma = 16; % SD in Gaussian description of peak signal

xx = 1:scanLength; % x-coordinate of scan

bckav = 25; % Average background signal

sigpk = 2; % Average signal peak above background

backSig = zeros(1,scanLength); % Initialized background signal

sig = zeros(1,scanLength); % Initialized peak signal

% Initialize plots of progress in fitting

sigPkArray = zeros(1,numScans);

sigPosArray = zeros(1,numScans);

sigSigmaArray = zeros(1,numScans);

bckArray = zeros(1,numScans);

for i = 1:numScans

subplot(4,5,[1,2,3,4,6,7,8,9,11,12,13,14,16,17,18,19])

newplot

hold off

% Simulated background with Poisson noise

backSig = backSig+poissrnd(bckav*ones(1,scanLength));

% Idealized peak signal

perfSig = sigpk*exp(-(xx - scanLength/2).^2/(2*sigma^2));

% Simulated peak signal with Poisson noise

sig = sig+poissrnd(perfSig);

% Add up simulated noisy signals and normalized by dividing by i

totSig = (sig + backSig)/i;

plot(xx,totSig,'color',philBlue,'LineWidth', 2);

set(gca,'linewidth',2,'fontsize',18);

ymin = 0.98*min(totSig);

ymax = 1.02*max(totSig);

ylim([ymin ymax]);

fun = @(x,xx)x(1)*exp(-(xx-x(2)).^2/(2*x(3)^2))+x(4);

x0 = [10,scanLength/2.2,sigma*2,10]; % Random starting guess

x = lsqcurvefit(fun,x0,xx,totSig);

hold on

fity = x(1)*exp(-(xx-x(2)).^2/(2*x(3)^2))+x(4);

plot(xx,fity,'color','red','LineWidth', 2);

set(gca,'linewidth',2,'fontsize',14);

str1 = 'Scan ';

str2 = num2str(i,'% 4.0f');

strTot = [str1,str2];

hTextYpos = ymin+(ymax-ymin)*0.88;

hText = text(0.88*scanLength, hTextYpos, strTot, 'FontSize',18);

sigPkArray(i) = x(1);

sigPosArray(i) = x(2);

sigSigmaArray(i) = x(3);

bckArray(i) = x(4);

%hold on

subplot(4,5,[5])

plot(sigPkArray(1:i),'color','red','LineWidth', 2);

set(gca,'linewidth',2,'fontsize',14);

xlim([1 numScans]);

ylim([sigpk*0.5 sigpk*1.5]);

hText = text(numScans*0.8, sigpk*1.25, 'a_0', 'FontSize',14);

subplot(4,5,[10])

plot(sigPosArray(1:i),'color','red','LineWidth', 2);

set(gca,'linewidth',2,'fontsize',14);

xlim([1 numScans]);

ylim([sigPos*0.98 sigPos*1.02]);

hText = text(numScans*0.8, sigPos*1.01, 'a_1', 'FontSize',14);

subplot(4,5,[15])

plot(sigSigmaArray(1:i),'color','red','LineWidth', 2);

set(gca,'linewidth',2,'fontsize',14);

xlim([1 numScans]);

ylim([sigma*0.5 sigma*1.5]);

hText = text(numScans*0.8, sigma*1.25, 'a_2', 'FontSize',14);

subplot(4,5,[20])

plot(bckArray(1:i),'color','red','LineWidth', 2);

set(gca,'linewidth',2,'fontsize',14);

xlim([1 numScans]);

ylim([bckav*0.98 bckav*1.02]);

hText = text(numScans*0.8, bckav*1.01, 'a_3', 'FontSize',14);

% Store the frame

frame = getframe(gcf);

writeVideo(vid,frame);

end

% Output the movie as an mpg file

close(vid);