diff --git a/OpenSIM manual.pdf b/OpenSIM manual.pdf
new file mode 100644
index 0000000..49b26d9
Binary files /dev/null and b/OpenSIM manual.pdf differ
diff --git a/SIMbasic/ApodizationFunction.m b/SIMbasic/ApodizationFunction.m
new file mode 100644
index 0000000..172f322
--- /dev/null
+++ b/SIMbasic/ApodizationFunction.m
@@ -0,0 +1,30 @@
+function [FsumA] = ApodizationFunction(Fsum,k2fa,k2fb,k2fc,Kotf,Index)
+
+t = size(Fsum,1);
+
+OTFo_mask = OTFmaskShifted(0.*k2fa,Kotf,t);
+OTFap_mask = OTFmaskShifted(k2fa,Kotf,t);
+OTFam_mask = OTFmaskShifted(-k2fa,Kotf,t);
+OTFbp_mask = OTFmaskShifted(k2fb,Kotf,t);
+OTFbm_mask = OTFmaskShifted(-k2fb,Kotf,t);
+OTFcp_mask = OTFmaskShifted(k2fc,Kotf,t);
+OTFcm_mask = OTFmaskShifted(-k2fc,Kotf,t);
+ApoMask = OTFo_mask.*OTFap_mask.*OTFam_mask.*OTFbp_mask.*OTFbm_mask.*OTFcp_mask.*OTFcm_mask;
+clear OTFoz OTFo_mask OTFap_mask OTFam_mask OTFbp_mask OTFbm_mask OTFcp_mask OTFcm_mask
+%{
+figure;
+imshow(ApoMask,[ ])
+%}
+
+DistApoMask = bwdist(ApoMask);
+maxApoMask = max(max(DistApoMask));
+ApoFunc = double(DistApoMask./maxApoMask).^Index;
+%{
+figure;
+imshow(DistApoMask,[ ])
+figure;
+mesh(double(ApoFunc))
+%}
+
+FsumA = Fsum.*ApoFunc;
+clear ApoMask DistApoMask ApoFunc
\ No newline at end of file
diff --git a/SIMbasic/ApproxFreqDuplex.m b/SIMbasic/ApproxFreqDuplex.m
new file mode 100644
index 0000000..b30640d
--- /dev/null
+++ b/SIMbasic/ApproxFreqDuplex.m
@@ -0,0 +1,30 @@
+function [maxK2,Ix,Iy] = ApproxFreqDuplex(FiSMap,Kotf)
+% AIM: approx. illumination frequency vector determination
+% INPUT VARIABLES
+% FiSMap: FT of raw SIM image
+% Kotf: OTF cut-off frequency
+% OUTPUT VARIABLES
+% maxK2: illumination frequency vector (approx)
+% Ix,Iy: coordinates of illumination frequency peaks
+
+FiSMap = abs(FiSMap);
+
+w = size(FiSMap,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+Ro = sqrt( (X-wo).^2 + (Y-wo).^2 );
+Z0 = Ro > round(0.5*Kotf);
+Z1 = X > wo;
+
+FiSMap = FiSMap.*Z0.*Z1;
+dumY = max( FiSMap,[],1 );
+[dummy Iy] = max(dumY);
+dumX = max( FiSMap,[],2 );
+[dummy Ix] = max(dumX);
+
+maxK2 = [Ix-(wo+1) Iy-(wo+1)];
+
+
diff --git a/SIMbasic/IlluminationFreqF.m b/SIMbasic/IlluminationFreqF.m
new file mode 100644
index 0000000..8dd5ec4
--- /dev/null
+++ b/SIMbasic/IlluminationFreqF.m
@@ -0,0 +1,38 @@
+function [k2fa] = IlluminationFreqF(S1aTnoisy,OTFo)
+% AIM: illumination frequency vector determination
+% INPUT VARIABLES
+% S1aTnoisy: raw SIM image
+% OTFo: system OTF
+% OUTPUT VARIABLE
+% k2fa: illumination frequency vector
+
+w = size(OTFo,1);
+wo = w/2;
+
+% computing PSFe for edge tapering SIM images
+PSFd = real(fftshift( ifft2(fftshift(OTFo.^10)) ));
+PSFd = PSFd/max(max(PSFd));
+PSFd = PSFd/sum(sum(PSFd));
+h = 30;
+PSFe = PSFd(wo-h+1:wo+h,wo-h+1:wo+h);
+
+% edge tapering raw SIM image
+S1aTnoisy_et = edgetaper(S1aTnoisy,PSFe);
+fS1aTnoisy_et = fftshift(fft2(S1aTnoisy_et));
+
+% OTF cut-off freq
+Kotf = OTFedgeF(OTFo);
+
+% Approx illumination frequency vector
+[k2fa,~,~] = ApproxFreqDuplex(fS1aTnoisy_et,Kotf);
+
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+% illumination frequency vector determination by optimizing
+% autocorrelation of fS1aTnoisy
+OPT = 1;
+PhaseKai2opt0 = @(k2fa0)PhaseKai2opt(k2fa0,fS1aTnoisy,OTFo,OPT);
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+k2fa0 = k2fa;
+[k2fa,fval] = fminsearch(PhaseKai2opt0,k2fa0,options);
+% k2a = sqrt(k2fa*k2fa')
diff --git a/SIMbasic/IlluminationPhaseF.m b/SIMbasic/IlluminationPhaseF.m
new file mode 100644
index 0000000..2898568
--- /dev/null
+++ b/SIMbasic/IlluminationPhaseF.m
@@ -0,0 +1,14 @@
+function [phaseA1] = IlluminationPhaseF(S1aTnoisy,k2fa)
+% AIM: illumination phase shift determination
+% INPUT VARIABLES
+% S1aTnoisy: raw SIM image
+% k2fa: illumination frequency vector
+% OUTPUT VARIABLE
+% phaseA1: illumination phase shift determined
+
+PatternPhaseOpt0 = @(phaseA0)PatternPhaseOpt(phaseA0,S1aTnoisy,k2fa);
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+phaseA0 = 0;
+[phaseA1,fval] = fminsearch(PatternPhaseOpt0,phaseA0,options);
+% phaseA1*180/pi
\ No newline at end of file
diff --git a/SIMbasic/MergingHeptaletsF.m b/SIMbasic/MergingHeptaletsF.m
new file mode 100644
index 0000000..19a98f6
--- /dev/null
+++ b/SIMbasic/MergingHeptaletsF.m
@@ -0,0 +1,68 @@
+function [Fsum,Fperi,Fcent] = MergingHeptaletsF(fDIo,fDIp,fDIm,fDBo,fDBp,fDBm,...
+ fDCo,fDCp,fDCm,Ma,Mb,Mc,npDIo,npDIp,npDIm,...
+ npDBo,npDBp,npDBm,npDCo,npDCp,npDCm,k2fa,k2fb,k2fc,OBJpara,OTFo)
+
+% AIM: To merge all 9 frequency components into one using generalized
+% Wiener Filter
+% INPUT VARIABLES
+% [fDIo fDIp fDIm
+% fDBo fDBp fDBm
+% fDCo fDCp fDCm]: nine frequency components
+% Ma,Mb,Mc: modulation factors for the three illumination orientations
+% [npDIo npDIp npDIm
+% npDBo npDBp npDBm
+% npDCo,npDCp,npDCm]: noise powers corresponding to nine frequency
+% components
+% k2fa,k2fb,k2fc: illumination frequency vectors for the three
+% illumination orientations
+% OBJpara: Object spectrum parameters
+% OTFo: system OTF
+% OUTPUT VARIABLES
+% Fsum: all nine frequency components merged into one using
+% generalised Wiener Filter
+% Fperi: six off-center frequency components merged into one using
+% generalised Wiener Filter
+% Fcent: averaged of the three central frequency components
+
+
+% obtain signal spectrums corresponding to central and off-center
+% frequency components
+[SIGao,SIGam2,SIGap2] = TripletSNR0(OBJpara,k2fa,OTFo,fDIm);
+[SIGbo,SIGbm2,SIGbp2] = TripletSNR0(OBJpara,k2fb,OTFo,fDBm);
+[SIGco,SIGcm2,SIGcp2] = TripletSNR0(OBJpara,k2fc,OTFo,fDCm);
+
+SIGap2 = Ma*SIGap2;
+SIGam2 = Ma*SIGam2;
+SIGbp2 = Mb*SIGbp2;
+SIGbm2 = Mb*SIGbm2;
+SIGcp2 = Mc*SIGcp2;
+SIGcm2 = Mc*SIGcm2;
+
+%% Generalized Wiener-Filter computation
+SNRao = SIGao.*conj(SIGao)./npDIo;
+SNRap = SIGap2.*conj(SIGap2)./npDIp;
+SNRam = SIGam2.*conj(SIGam2)./npDIm;
+
+SNRbo = SIGbo.*conj(SIGbo)./npDBo;
+SNRbp = SIGbp2.*conj(SIGbp2)./npDBp;
+SNRbm = SIGbm2.*conj(SIGbm2)./npDBm;
+
+SNRco = SIGco.*conj(SIGco)./npDCo;
+SNRcp = SIGcp2.*conj(SIGcp2)./npDCp;
+SNRcm = SIGcm2.*conj(SIGcm2)./npDCm;
+
+ComDeno = 0.01 + ( SNRao + SNRap + SNRam + SNRbo + SNRbp + SNRbm + SNRco + SNRcp + SNRcm );
+Fsum = fDIo.*SNRao + fDIp.*SNRap + fDIm.*SNRam...
+ + fDBo.*SNRbo + fDBp.*SNRbp + fDBm.*SNRbm...
+ + fDCo.*SNRco + fDCp.*SNRcp + fDCm.*SNRcm;
+Fsum = Fsum./ComDeno;
+
+ComPeri = 0.01 + ( SNRap + SNRam + SNRbp + SNRbm + SNRcp + SNRcm );
+Fperi = fDIp.*SNRap + fDIm.*SNRam + fDBp.*SNRbp + fDBm.*SNRbm...
+ + fDCp.*SNRcp + fDCm.*SNRcm;
+Fperi = Fperi./ComPeri;
+
+% averaged central frequency component
+Fcent = (fDIo+fDBo+fDCo)/3;
+
+
diff --git a/SIMbasic/ModulationFactor.m b/SIMbasic/ModulationFactor.m
new file mode 100644
index 0000000..2ac2dda
--- /dev/null
+++ b/SIMbasic/ModulationFactor.m
@@ -0,0 +1,86 @@
+function Mm = ModulationFactor(fDp,k1a,OBJparaA,OTFo)
+% AIM: Determination of modulation factor
+% INPUT VARIABLES
+% fDp: off-center frequency component
+% k1a: illumination frequency vector
+% OBJparaA: Object power parameters
+% OTFo: system OTF
+% OUTPUT VARIABLE
+% Mm: modulation factor
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% magnitude of illumination vector
+k2 = sqrt(k1a*k1a');
+
+% vector along illumination direction
+kv = k1a(2) + 1i*k1a(1);
+Rp = abs(Cv+kv);
+
+% Object power parameters
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+
+% Object spectrum
+OBJp = Aobj*(Rp+0).^Bobj;
+
+% illumination vector rounded to nearest pixel
+k3 = -round(k1a);
+
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+
+% signal spectrum
+SIGap = OBJp.*OTFo;
+
+% figure;
+% mesh(log(abs(OBJp)))
+% figure;
+% mesh(log(abs(fDp)))
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = fDp.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Noise free object power computation
+Fpower = fDp.*conj(fDp) - NoisePower;
+fDp = sqrt(abs(Fpower));
+
+% frequency range over which signal power matching is done to estimate
+% modulation factor
+Zmask = (Ro > 0.2*k2).*(Ro < 0.8*k2).*(Rp > 0.2*k2);
+
+% least square approximation for modulation factor
+Mm = sum(sum(SIGap.*abs(fDp).*Zmask));
+Mm = Mm./sum(sum(SIGap.^2.*Zmask))
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(fDp(wo+pp,:))),'o-','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(SIGap(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,log(abs(Mm*SIGap(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+grid on
+box on
+kk
+%}
+
+
+
diff --git a/SIMbasic/OBJparaOpt.m b/SIMbasic/OBJparaOpt.m
new file mode 100644
index 0000000..8acbbbb
--- /dev/null
+++ b/SIMbasic/OBJparaOpt.m
@@ -0,0 +1,47 @@
+function Esum = OBJparaOpt(OBJpara0,fDIoTnoisy,OTFo)
+% AIM: Determined Sum of Squared Errors (SSE) between `actual signal power' and
+% `approximated signal power'
+% INPUT VARIABLES
+% OBJpara0: [Aobj Bobj], Object power parameters
+% fDIoTnoisy: FT of central frequency component
+% OTFo: system OTF
+% OUTPUT VARIABLE
+% Esum: SSE between `actual signal spectrum' and `approximated signal spectrum'
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+Ro(wo+1,wo+1) = 1; % to avoid nan
+
+% approximated signal power calculation
+Aobj = OBJpara0(1);
+Bobj = OBJpara0(2);
+OBJpower = Aobj*(Ro.^Bobj);
+SIGpower = OBJpower.*OTFo;
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% range of frequency over which SSE is computed
+Zloop = (Ro<0.75*Kotf).*(Ro>0.25*Kotf);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = fDIoTnoisy.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Noise free object power computation
+Fpower = fDIoTnoisy.*conj(fDIoTnoisy) - NoisePower;
+fDIoTnoisy = sqrt(abs(Fpower));
+
+% SSE computation
+Error = fDIoTnoisy - SIGpower;
+Esum = sum(sum((Error.^2./Ro).*Zloop));
+
diff --git a/SIMbasic/OBJpowerPara.m b/SIMbasic/OBJpowerPara.m
new file mode 100644
index 0000000..9ef86ce
--- /dev/null
+++ b/SIMbasic/OBJpowerPara.m
@@ -0,0 +1,48 @@
+function OBJparaA = OBJpowerPara(fDIoTnoisy,OTFo)
+% AIM: determination of object power parameters Aobj and Bobj
+% INPUT VARIABLES
+% fDIoTnoisy: FT of central frequency component
+% OTFo: system OTFo
+% OUTPUT VARIABLES
+% OBJparaA: [Aobj Bobj]
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+%% object power parameters through optimization
+OBJparaOpt0 = @(OBJpara0)OBJparaOpt(OBJpara0,fDIoTnoisy,OTFo);
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+
+% obtaining crude initial guesses for Aobj and Bobj
+Zm = (Ro>0.3*Kotf).*(Ro<0.4*Kotf);
+Aobj = sum(sum(abs(fDIoTnoisy.*Zm)))./sum(sum(Zm));
+Bobj = -0.5;
+OBJpara0 = [Aobj Bobj];
+
+% optimization step
+[OBJparaA,fval] = fminsearch(OBJparaOpt0,OBJpara0,options);
+
+%% plot for cross-checking the result
+%{
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+OBJo = Aobj*(Ro.^Bobj);
+SIGp = OBJo.*OTFo;
+pp = 3;
+figure;
+hold on
+plot([0:w-1]-wo,log(abs(fDIoTnoisy(wo+pp,:))),'k--','LineWidth',3,'MarkerSize',6)
+plot([0:w-1]-wo,log(abs(SIGp(wo+pp,:))),'mo-','LineWidth',2,'MarkerSize',6)
+legend('acutal signal power','avg. signal power')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/SIMbasic/OTFdoubling.m b/SIMbasic/OTFdoubling.m
new file mode 100644
index 0000000..58aaa6e
--- /dev/null
+++ b/SIMbasic/OTFdoubling.m
@@ -0,0 +1,17 @@
+function OTF2 = OTFdoubling(OTFo,DoubleMatSize)
+
+%% embeds OTFo in doubled range of frequencies
+
+w = size(OTFo,1);
+wo = w/2;
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ OTF2 = zeros(2*w,2*w);
+ OTF2(wo+1:w+wo,wo+1:w+wo) = OTFo;
+ clear OTFo
+else
+ OTF2 = OTFo;
+ clear OTFo
+end
\ No newline at end of file
diff --git a/SIMbasic/OTFedgeF.m b/SIMbasic/OTFedgeF.m
new file mode 100644
index 0000000..e4d6a7f
--- /dev/null
+++ b/SIMbasic/OTFedgeF.m
@@ -0,0 +1,13 @@
+function Kotf = OTFedgeF(OTFo)
+
+w = size(OTFo,1);
+wo = w/2;
+
+OTF1 = OTFo(wo+1,:);
+OTFmax = max(max(abs(OTFo)));
+OTFtruncate = 0.01;
+i = 1;
+while ( abs(OTF1(1,i))<OTFtruncate*OTFmax )
+ Kotf = wo+1-i;
+ i = i + 1;
+end
\ No newline at end of file
diff --git a/SIMbasic/OTFmaskShifted.m b/SIMbasic/OTFmaskShifted.m
new file mode 100644
index 0000000..a60ac98
--- /dev/null
+++ b/SIMbasic/OTFmaskShifted.m
@@ -0,0 +1,13 @@
+function [OTFap_mask] = OTFmaskShifted(k2fa,Kotf,w)
+% OTFap_mask = complement of shifted OTF
+
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+Iy = k2fa(1,1) + (wo+1);
+Ix = k2fa(1,2) + (wo+1);
+
+Ro = sqrt( (X+1-Ix).^2 + (Y+1-Iy).^2 );
+OTFap_mask = Ro > Kotf;
\ No newline at end of file
diff --git a/SIMbasic/PCMfilteringF.m b/SIMbasic/PCMfilteringF.m
new file mode 100644
index 0000000..3a04408
--- /dev/null
+++ b/SIMbasic/PCMfilteringF.m
@@ -0,0 +1,143 @@
+function [fDof,fDp2,fDm2,npDo,npDp,npDm,Mm,DoubleMatSize]...
+ = PCMfilteringF(fDo,fDp,fDm,OTFo,OBJparaA,kA)
+
+% AIM: obtaining Wiener Filtered estimates of noisy frequency components
+% INPUT VARIABLES
+% fDo,fDp,fDm: noisy estimates of separated frequency components
+% OTFo: system OTF
+% OBJparaA: object power parameters
+% kA: illumination vector
+% OUTPUT VARIABLES
+% fDof,fDp2,fDm2: Wiener Filtered estimates of fDo,fDp,fDm
+% npDo,npDp,npDm: avg. noise power in fDo,fDp,fDm
+% Mm: illumination modulation factor
+% DoubleMatSize: parameter for doubling FT size if necessary
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% object power parameters
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+
+% Wiener Filtering central frequency component
+SFo = 1;
+co = 1.0;
+[fDof,npDo] = WoFilterCenter(fDo,OTFo,co,OBJparaA,SFo);
+
+% modulation factor determination
+Mm = ModulationFactor(fDp,kA,OBJparaA,OTFo);
+
+%% Duplex power (default)
+kv = kA(2) + 1i*kA(1); % vector along illumination direction
+Rp = abs(Cv-kv);
+Rm = abs(Cv+kv);
+OBJp = Aobj*(Rp.^Bobj);
+OBJm = Aobj*(Rm.^Bobj);
+k3 = round(kA);
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+OBJm(wo+1-k3(1),wo+1-k3(2)) = 0.25*OBJm(wo+2-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+2-k3(2))...
+ + 0.25*OBJm(wo+0-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+0-k3(2));
+
+% Filtering side lobes (off-center frequency components)
+SFo = Mm;
+[fDpf,npDp] = WoFilterSideLobe(fDp,OTFo,co,OBJm,SFo);
+[fDmf,npDm] = WoFilterSideLobe(fDm,OTFo,co,OBJp,SFo);
+
+%% doubling Fourier domain size if necessary
+DoubleMatSize = 0;
+if ( 2*Kotf > wo )
+ DoubleMatSize = 1; % 1 for doubling fourier domain size, 0 for keeping it unchanged
+end
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ to = t/2;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ [U,V] = meshgrid(u,v);
+ fDoTemp = zeros(2*w,2*w);
+ fDpTemp = zeros(2*w,2*w);
+ fDmTemp = zeros(2*w,2*w);
+ OTFtemp = zeros(2*w,2*w);
+ fDoTemp(wo+1:w+wo,wo+1:w+wo) = fDof;
+ fDpTemp(wo+1:w+wo,wo+1:w+wo) = fDpf;
+ fDmTemp(wo+1:w+wo,wo+1:w+wo) = fDmf;
+ OTFtemp(wo+1:w+wo,wo+1:w+wo) = OTFo;
+ clear fDof fDpf fDmf OTFo w wo x y X Y
+ fDof = fDoTemp;
+ fDpf = fDpTemp;
+ fDmf = fDmTemp;
+ OTFo = OTFtemp;
+ clear fDoTemp fDpTemp fDmTemp OTFtemp
+else
+ t = w;
+ to = t/2;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ [U,V] = meshgrid(u,v);
+end
+
+% Shifting the off-center frequency components to their correct location
+fDp1 = fft2(ifft2(fDpf).*exp( +1i.*2*pi*(kA(2)/t.*(U-to) + kA(1)/t.*(V-to)) ));
+fDm1 = fft2(ifft2(fDmf).*exp( -1i.*2*pi*(kA(2)/t.*(U-to) + kA(1)/t.*(V-to)) ));
+
+%% Shift induced phase error correction
+Cv = (U-to) + 1i*(V-to);
+Ro = abs(Cv);
+Rp = abs(Cv-kv);
+k2 = sqrt(kA*kA');
+
+% frequency range over which corrective phase is determined
+Zmask = (Ro < 0.8*k2).*(Rp < 0.8*k2);
+
+% corrective phase
+Angle0 = angle( sum(sum( fDof.*conj(fDp1).*Zmask )) );
+
+% phase correction
+fDp2 = exp(+1i*Angle0).*fDp1;
+fDm2 = exp(-1i*Angle0).*fDm1;
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot(u-to,angle(fDof(to+pp,:)).*180/pi,'o-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(fDp1(to+pp,:)).*180/pi,'ro-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(fDp2(to+pp,:)).*180/pi,'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+figure;
+hold on
+plot3(u-to,real(fDof(to+1,:)),imag(fDof(to+1,:)),'o-','LineWidth',2,'MarkerSize',8)
+plot3(u-to,real(fDp1(to+1,:)),imag(fDp1(to+1,:)),'ro-','LineWidth',2,'MarkerSize',8)
+plot3(u-to,real(fDp2(to+1,:)),imag(fDp2(to+1,:)),'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+
+ff1 = fDof.*conj(fDp1);
+ff2 = fDof.*conj(fDp2);
+figure;
+hold on
+plot(u-to,angle(ff1(to+pp,:)).*180/pi,'ro-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(ff2(to+pp,:)).*180/pi,'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+%}
+
+
+
+
diff --git a/SIMbasic/PCMseparateF.m b/SIMbasic/PCMseparateF.m
new file mode 100644
index 0000000..d6b12d8
--- /dev/null
+++ b/SIMbasic/PCMseparateF.m
@@ -0,0 +1,64 @@
+function [fDo,fDp,fDm,kA]...
+ = PCMseparateF(S1aTnoisy,S2aTnoisy,S3aTnoisy,OTFo)
+% AIM: obtaining the noisy estimates of three frequency components
+% INPUT VARIABLES
+% S1aTnoisy,S2aTnoisy,S3aTnoisy: 3 raw SIM images with identical
+% illumination pattern orientation
+% but different phase shifts
+% OTFo: system OTF
+% OUTPUT VARIABLES
+% fDo,fDp,fDm: noisy estimates of separated frequency components
+% kA: (averaged) illumination frequency vector
+
+
+w = size(S1aTnoisy,1);
+wo = w/2;
+
+%% Determination of illumination frequency vectors
+[k1a] = IlluminationFreqF(S1aTnoisy,OTFo);
+[k2a] = IlluminationFreqF(S2aTnoisy,OTFo);
+[k3a] = IlluminationFreqF(S3aTnoisy,OTFo);
+
+% mean illumination frequency vector
+kA = (k1a + k2a + k3a)/3;
+
+%% determination of illumination phase shifts
+[phase1A] = IlluminationPhaseF(S1aTnoisy,kA);
+[phase2A] = IlluminationPhaseF(S2aTnoisy,kA);
+[phase3A] = IlluminationPhaseF(S3aTnoisy,kA);
+
+% for display in command window
+phaseA = [phase1A; phase2A; phase3A];
+phaseA*180/pi
+
+% computing PSFe for edge tapering SIM images
+PSFd = real(fftshift( ifft2(fftshift(OTFo.^3)) ));
+PSFd = PSFd/max(max(PSFd));
+PSFd = PSFd/sum(sum(PSFd));
+h = 30;
+PSFe = PSFd(wo-h+1:wo+h,wo-h+1:wo+h);
+
+% edge tapering raw SIM images
+S1aTnoisy = edgetaper(S1aTnoisy,PSFe);
+S2aTnoisy = edgetaper(S2aTnoisy,PSFe);
+S3aTnoisy = edgetaper(S3aTnoisy,PSFe);
+
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+fS2aTnoisy = fftshift(fft2(S2aTnoisy));
+fS3aTnoisy = fftshift(fft2(S3aTnoisy));
+
+%% Separating the three frequency components
+phaseShift0 = phase1A;
+phaseShift = [phase2A; phase3A];
+[fDo,fDp,fDm] = SeparatedComponents2D(...
+ phaseShift,phaseShift0,fS1aTnoisy,fS2aTnoisy,fS3aTnoisy);
+
+%% for visual verification
+%{
+figure;
+mesh(log(abs(fDo)))
+figure;
+mesh(log(abs(fDp)))
+figure;
+mesh(log(abs(fDm)))
+%}
diff --git a/SIMbasic/PatternPhaseOpt.m b/SIMbasic/PatternPhaseOpt.m
new file mode 100644
index 0000000..022539a
--- /dev/null
+++ b/SIMbasic/PatternPhaseOpt.m
@@ -0,0 +1,12 @@
+function CCop = PatternPhaseOpt(phaseA,S1aTnoisy,k2fa)
+
+w = size(S1aTnoisy,1);
+wo = w/2;
+
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+sAo = cos( 2*pi*(k2fa(2).*(X-wo)+k2fa(1).*(Y-wo))./w + phaseA );
+S1aTnoisy = S1aTnoisy - mean2(S1aTnoisy);
+CCop = -sum(sum(S1aTnoisy.*sAo));
\ No newline at end of file
diff --git a/SIMbasic/PhaseKai2opt.m b/SIMbasic/PhaseKai2opt.m
new file mode 100644
index 0000000..bb7ec1f
--- /dev/null
+++ b/SIMbasic/PhaseKai2opt.m
@@ -0,0 +1,46 @@
+function [CCop] = PhaseKai2opt(k2fa,fS1aTnoisy,OTFo,OPT)
+% Aim: Compute autocorrelation of FT of raw SIM images
+% k2fa: illumination frequency vector
+% fS1aTnoisy: FT of raw SIM image
+% OTFo: system OTF
+% OPT: acronym for `OPTIMIZE'; to be set to 1 when this function is used
+% for optimization, or else to 0
+% CCop: autocorrelation of fS1aTnoisy
+
+w = size(fS1aTnoisy,1);
+wo = w/2;
+
+fS1aTnoisy = fS1aTnoisy.*(1-1*OTFo.^10);
+fS1aT = fS1aTnoisy.*conj(OTFo);
+
+Kotf = OTFedgeF(OTFo);
+DoubleMatSize = 0;
+if ( 2*Kotf > wo )
+ DoubleMatSize = 1; % 1 for doubling fourier domain size, 0 for keeping it unchanged
+end
+
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ fS1aT_temp = zeros(t,t);
+ fS1aT_temp(wo+1:w+wo,wo+1:w+wo) = fS1aT;
+ clear fS1aT
+ fS1aT = fS1aT_temp;
+ clear fS1aT_temp
+else
+ t = w;
+end
+to = t/2;
+u = linspace(0,t-1,t);
+v = linspace(0,t-1,t);
+[U,V] = meshgrid(u,v);
+S1aT = exp( -1i.*2*pi*( k2fa(2)/t.*(U-to)+k2fa(1)/t.*(V-to) ) ).*ifft2(fS1aT);
+fS1aT0 = fft2( S1aT );
+
+mA = sum(sum( fS1aT.*conj(fS1aT0) ));
+mA = mA./sum(sum( fS1aT0.*conj(fS1aT0) ));
+
+if (OPT > 0)
+ CCop = -abs(mA);
+else
+ CCop = mA;
+end
\ No newline at end of file
diff --git a/SIMbasic/PsfOtf.m b/SIMbasic/PsfOtf.m
new file mode 100644
index 0000000..0120d1c
--- /dev/null
+++ b/SIMbasic/PsfOtf.m
@@ -0,0 +1,34 @@
+function [yy0,OTF2dc] = PsfOtf(w,scale)
+
+% AIM: To generate PSF and OTF using Bessel function
+% INPUT VARIABLES
+% w: image size
+% scale: a parameter used to adjust PSF/OTF width
+% OUTPUT VRAIBLES
+% yyo: system PSF
+% OTF2dc: system OTF
+
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+%% Generation of the PSF with Besselj.
+% scale = 0.3;
+R=sqrt(min(X,abs(X-w)).^2+min(Y,abs(Y-w)).^2);
+yy=abs(2*besselj(1,scale*R+eps,1)./(scale*R+eps)).^2; %0.5 is introduced to make PSF wider
+yy0 = fftshift(yy);
+% figure;
+% mesh(X,Y,yy0)
+% figure;
+% imshow(yy0,[])
+% axis square
+
+%Generate 2D OTF.
+OTF2d=fft2(yy);
+OTF2dmax = max(max(abs(OTF2d)));
+OTF2d = OTF2d./OTF2dmax;
+OTF2dc = abs(fftshift(OTF2d));
+% figure;
+% mesh(abs(OTF2dc))
+% title('OTF of the system')
+% axis square
\ No newline at end of file
diff --git a/SIMbasic/SIMimagesF.m b/SIMbasic/SIMimagesF.m
new file mode 100644
index 0000000..6b93dba
--- /dev/null
+++ b/SIMbasic/SIMimagesF.m
@@ -0,0 +1,282 @@
+function [S1aTnoisy S2aTnoisy S3aTnoisy ...
+ S1bTnoisy S2bTnoisy S3bTnoisy ...
+ S1cTnoisy S2cTnoisy S3cTnoisy ...
+ DIoTnoisy DIoT] = SIMimagesF(k2,...
+ DIo,PSFo,OTFo,ModFac,NoiseLevel,UsePSF)
+
+% AIM: to generate raw sim images
+% INPUT VARIABLES
+% k2: illumination frequency
+% DIo: specimen image
+% PSFo: system PSF
+% OTFo: system OTF
+% UsePSF: 1 (to blur SIM images by convloving with PSF)
+% 0 (to blur SIM images by truncating its fourier content beyond OTF)
+% NoiseLevel: percentage noise level for generating gaussian noise
+% OUTPUT VARIABLES
+% [S1aTnoisy S2aTnoisy S3aTnoisy
+% S1bTnoisy S2bTnoisy S3bTnoisy
+% S1cTnoisy S2cTnoisy S3cTnoisy]: nine raw sim images
+% DIoTnoisy: noisy wide field image
+% DIoT: noise-free wide field image
+
+w = size(DIo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+%% illunination phase shifts along the three directions
+p0Ao = 0*pi/3;
+p0Ap = 2*pi/3;
+p0Am = 4*pi/3;
+p0Bo = 0*pi/3;
+p0Bp = 2*pi/3;
+p0Bm = 4*pi/3;
+p0Co = 0*pi/3;
+p0Cp = 2*pi/3;
+p0Cm = 4*pi/3;
+
+%% Illuminating patterns
+alpha = 0*pi/6;
+% orientation direction of illumination patterns
+thetaA = 0*pi/3 + alpha;
+thetaB = 1*pi/3 + alpha;
+thetaC = 2*pi/3 + alpha;
+% illumination frequency vectors
+k2a = (k2/w).*[cos(thetaA) sin(thetaA)];
+k2b = (k2/w).*[cos(thetaB) sin(thetaB)];
+k2c = (k2/w).*[cos(thetaC) sin(thetaC)];
+% -------------------------------------------------------
+% mean illumination intensity
+mA = 0.5;
+mB = 0.5;
+mC = 0.5;
+% amplitude of illumination intensity above mean
+aA = 0.5*ModFac;
+aB = 0.5*ModFac;
+aC = 0.5*ModFac;
+
+% random phase shift errors
+NN = 1*(0.5-rand(9,1))*pi/18;
+
+% illunination phase shifts with random errors
+psAo = p0Ao + NN(1,1);
+psAp = p0Ap + NN(2,1);
+psAm = p0Am + NN(3,1);
+psBo = p0Bo + NN(4,1);
+psBp = p0Bp + NN(5,1);
+psBm = p0Bm + NN(6,1);
+psCo = p0Co + NN(7,1);
+psCp = p0Cp + NN(8,1);
+psCm = p0Cm + NN(9,1);
+
+% illunination patterns
+sAo = mA + aA*cos(2*pi*(k2a(1,1).*(X-wo)+k2a(1,2).*(Y-wo))+psAo); % illuminated signal (0 phase)
+sAp = mA + aA*cos(2*pi*(k2a(1,1).*(X-wo)+k2a(1,2).*(Y-wo))+psAp); % illuminated signal (+ phase)
+sAm = mA + aA*cos(2*pi*(k2a(1,1).*(X-wo)+k2a(1,2).*(Y-wo))+psAm); % illuminated signal (- phase)
+sBo = mB + aB*cos(2*pi*(k2b(1,1).*(X-wo)+k2b(1,2).*(Y-wo))+psBo); % illuminated signal (0 phase)
+sBp = mB + aB*cos(2*pi*(k2b(1,1).*(X-wo)+k2b(1,2).*(Y-wo))+psBp); % illuminated signal (+ phase)
+sBm = mB + aB*cos(2*pi*(k2b(1,1).*(X-wo)+k2b(1,2).*(Y-wo))+psBm); % illuminated signal (- phase)
+sCo = mC + aC*cos(2*pi*(k2c(1,1).*(X-wo)+k2c(1,2).*(Y-wo))+psCo); % illuminated signal (0 phase)
+sCp = mC + aC*cos(2*pi*(k2c(1,1).*(X-wo)+k2c(1,2).*(Y-wo))+psCp); % illuminated signal (+ phase)
+sCm = mC + aC*cos(2*pi*(k2c(1,1).*(X-wo)+k2c(1,2).*(Y-wo))+psCm); % illuminated signal (- phase)
+
+%{
+fftemp = fftshift(fft2(sAo));
+figure;
+hold on
+plot(x-wo,abs(OTFo(wo+1,:)),'--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,0.5*abs(fftemp(wo+1,:))./max(max(abs(fftemp))),'r--','LineWidth',2,'MarkerSize',6)
+legend('OTFo','illumination spectrum')
+grid on
+box on
+%}
+
+% figure;
+% imshow(sAo,[ ])
+% kk
+
+%% superposed Objects
+s1a = DIo.*sAo; % superposed signal (0 phase)
+s2a = DIo.*sAp; % superposed signal (+ phase)
+s3a = DIo.*sAm; % superposed signal (- phase)
+s1b = DIo.*sBo;
+s2b = DIo.*sBp;
+s3b = DIo.*sBm;
+s1c = DIo.*sCo;
+s2c = DIo.*sCp;
+s3c = DIo.*sCm;
+%{
+figure;
+subplot(3,3,1)
+imshow(s1a,[ ])
+title('s1a')
+subplot(3,3,4)
+imshow(s2a,[ ])
+title('s2a')
+subplot(3,3,7)
+imshow(s3a,[ ])
+title('s3a')
+subplot(3,3,2)
+imshow(s1b,[ ])
+title('s1b')
+subplot(3,3,5)
+imshow(s2b,[ ])
+title('s2b')
+subplot(3,3,8)
+imshow(s3b,[ ])
+title('s3b')
+subplot(3,3,3)
+imshow(s1c,[ ])
+title('s1c')
+subplot(3,3,6)
+imshow(s2c,[ ])
+title('s2c')
+subplot(3,3,9)
+imshow(s3c,[ ])
+title('s3c')
+%}
+
+
+%% superposed (noise-free) Images
+PSFsum = sum(sum(PSFo));
+if ( UsePSF == 1 )
+ DIoT = conv2(DIo,PSFo,'same')./PSFsum;
+ S1aT = conv2(s1a,PSFo,'same')./PSFsum;
+ S2aT = conv2(s2a,PSFo,'same')./PSFsum;
+ S3aT = conv2(s3a,PSFo,'same')./PSFsum;
+ S1bT = conv2(s1b,PSFo,'same')./PSFsum;
+ S2bT = conv2(s2b,PSFo,'same')./PSFsum;
+ S3bT = conv2(s3b,PSFo,'same')./PSFsum;
+ S1cT = conv2(s1c,PSFo,'same')./PSFsum;
+ S2cT = conv2(s2c,PSFo,'same')./PSFsum;
+ S3cT = conv2(s3c,PSFo,'same')./PSFsum;
+else
+ DIoT = ifft2( fft2(DIo).*fftshift(OTFo) );
+ S1aT = ifft2( fft2(s1a).*fftshift(OTFo) );
+ S2aT = ifft2( fft2(s2a).*fftshift(OTFo) );
+ S3aT = ifft2( fft2(s3a).*fftshift(OTFo) );
+ S1bT = ifft2( fft2(s1b).*fftshift(OTFo) );
+ S2bT = ifft2( fft2(s2b).*fftshift(OTFo) );
+ S3bT = ifft2( fft2(s3b).*fftshift(OTFo) );
+ S1cT = ifft2( fft2(s1c).*fftshift(OTFo) );
+ S2cT = ifft2( fft2(s2c).*fftshift(OTFo) );
+ S3cT = ifft2( fft2(s3c).*fftshift(OTFo) );
+
+ DIoT = real(DIoT);
+ S1aT = real(S1aT);
+ S2aT = real(S2aT);
+ S3aT = real(S3aT);
+ S1bT = real(S1bT);
+ S2bT = real(S2bT);
+ S3bT = real(S3bT);
+ S1cT = real(S1cT);
+ S2cT = real(S2cT);
+ S3cT = real(S3cT);
+end
+%{
+figure;
+subplot(3,3,1)
+imshow(S1aT,[ ])
+title('S1aT')
+subplot(3,3,4)
+imshow(S2aT,[ ])
+title('S2aT')
+subplot(3,3,7)
+imshow(S3aT,[ ])
+title('S3aT')
+subplot(3,3,2)
+imshow(S1bT,[ ])
+title('S1bT')
+subplot(3,3,5)
+imshow(S2bT,[ ])
+title('S2bT')
+subplot(3,3,8)
+imshow(S3bT,[ ])
+title('S3bT')
+subplot(3,3,3)
+imshow(S1cT,[ ])
+title('S1cT')
+subplot(3,3,6)
+imshow(S2cT,[ ])
+title('S2cT')
+subplot(3,3,9)
+imshow(S3cT,[ ])
+title('S3cT')
+%}
+
+
+%% Gaussian noise generation
+aNoise = NoiseLevel/100; % corresponds to 10% noise
+% SNR = 1/aNoise
+% SNRdb = 20*log10(1/aNoise)
+nDIoT = random('norm', 0, aNoise*std2(DIoT), w , w);
+nS1aT = random('norm', 0, aNoise*std2(S1aT), w , w);
+nS2aT = random('norm', 0, aNoise*std2(S2aT), w , w);
+nS3aT = random('norm', 0, aNoise*std2(S3aT), w , w);
+nS1bT = random('norm', 0, aNoise*std2(S2bT), w , w);
+nS2bT = random('norm', 0, aNoise*std2(S2bT), w , w);
+nS3bT = random('norm', 0, aNoise*std2(S2bT), w , w);
+nS1cT = random('norm', 0, aNoise*std2(S3cT), w , w);
+nS2cT = random('norm', 0, aNoise*std2(S3cT), w , w);
+nS3cT = random('norm', 0, aNoise*std2(S3cT), w , w);
+
+%% noise added raw SIM images
+NoiseFrac = 1; %may be set to 0 to avoid noise addition
+DIoTnoisy = DIoT + NoiseFrac*nDIoT;
+S1aTnoisy = S1aT + NoiseFrac*nS1aT;
+S2aTnoisy = S2aT + NoiseFrac*nS2aT;
+S3aTnoisy = S3aT + NoiseFrac*nS3aT;
+S1bTnoisy = S1bT + NoiseFrac*nS1bT;
+S2bTnoisy = S2bT + NoiseFrac*nS2bT;
+S3bTnoisy = S3bT + NoiseFrac*nS3bT;
+S1cTnoisy = S1cT + NoiseFrac*nS1cT;
+S2cTnoisy = S2cT + NoiseFrac*nS2cT;
+S3cTnoisy = S3cT + NoiseFrac*nS3cT;
+%{
+figure;
+subplot(3,3,1)
+imshow(S1aTnoisy,[ ])
+title('S1aTnoisy')
+subplot(3,3,4)
+imshow(S2aTnoisy,[ ])
+title('S2aTnoisy')
+subplot(3,3,7)
+imshow(S3aTnoisy,[ ])
+title('S3aTnoisy')
+subplot(3,3,2)
+imshow(S1bTnoisy,[ ])
+title('S1bTnoisy')
+subplot(3,3,5)
+imshow(S2bTnoisy,[ ])
+title('S2bTnoisy')
+subplot(3,3,8)
+imshow(S3bTnoisy,[ ])
+title('S3bTnoisy')
+subplot(3,3,3)
+imshow(S1cTnoisy,[ ])
+title('S1cTnoisy')
+subplot(3,3,6)
+imshow(S2cTnoisy,[ ])
+title('S2cTnoisy')
+subplot(3,3,9)
+imshow(S3cTnoisy,[ ])
+title('S3cTnoisy')
+%}
+
+%{
+figure;
+imshow(S1aT,[ ])
+title('S1aT')
+figure;
+imshow(S1aTnoisy,[ ])
+title('S1aTnoisy')
+
+figure;
+hold on
+plot(x-wo,S1aT(wo+1,:),'o-','LineWidth',2,'MarkerSize',6)
+plot(x-wo,S1aTnoisy(wo+1,:),'r*--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/SIMbasic/SIMo.m b/SIMbasic/SIMo.m
new file mode 100644
index 0000000..765b5e6
--- /dev/null
+++ b/SIMbasic/SIMo.m
@@ -0,0 +1,206 @@
+clear all
+close all
+clc
+
+w = 512;
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+%% Generation of the PSF with Besselj.
+scale = 0.63; % used to adjust PSF/OTF width
+[PSFo,OTFo] = PsfOtf(w,scale);
+%{
+figure;
+mesh(PSFo)
+title('system PSF')
+
+figure;
+mesh(OTFo)
+title('system OTF')
+%}
+
+%% Reading input file
+Io1 = imread('testpat.tiff');
+Io = Io1(257:768,257:768); % selecting central 512x512 of image
+DIo = double(Io);
+figure;
+imshow(DIo,[])
+
+
+%% Generating raw SIM Images
+k2 = 75.23; % illumination freq
+ModFac = 0.8; % modulation factor
+NoiseLevel = 10; % in percentage
+UsePSF = 0; % 1(to blur using PSF), 0(to blur using OTF)
+[S1aTnoisy S2aTnoisy S3aTnoisy ...
+ S1bTnoisy S2bTnoisy S3bTnoisy ...
+ S1cTnoisy S2cTnoisy S3cTnoisy ...
+ DIoTnoisy DIoT] = SIMimagesF(k2,DIo,PSFo,OTFo,ModFac,NoiseLevel,UsePSF);
+%{
+figure;
+imshow(S1aTnoisy,[])
+title('Raw SIM image')
+figure;
+imshow(S1bTnoisy,[])
+title('Raw SIM image')
+figure;
+imshow(S1cTnoisy,[])
+title('Raw SIM image')
+figure;
+imshow(DIoTnoisy,[])
+title('Noisy Wide-Field image')
+figure;
+imshow(DIoT,[])
+title('Noise-free Wide-Field image')
+break
+%}
+
+%% obtaining the noisy estimates of three frequency components
+[fAo,fAp,fAm,kA]...
+ = PCMseparateF(S1aTnoisy,S2aTnoisy,S3aTnoisy,OTFo);
+% break
+[fBo,fBp,fBm,kB]...
+ = PCMseparateF(S1bTnoisy,S2bTnoisy,S3bTnoisy,OTFo);
+[fCo,fCp,fCm,kC]...
+ = PCMseparateF(S1cTnoisy,S2cTnoisy,S3cTnoisy,OTFo);
+
+% averaging the central frequency components
+fCent = (fAo + fBo + fCo)/3;
+
+% Object power parameters determination
+OBJparaA = OBJpowerPara(fCent,OTFo);
+
+
+%% Wiener Filtering the noisy frequency components
+[fAof,fApf,fAmf,Nao,Nap,Nam,Ma,DoubleMatSize]...
+ = PCMfilteringF(fAo,fAp,fAm,OTFo,OBJparaA,kA);
+[fBof,fBpf,fBmf,Nbo,Nbp,Nbm,Mb,~]...
+ = PCMfilteringF(fBo,fBp,fBm,OTFo,OBJparaA,kB);
+[fCof,fCpf,fCmf,Nco,Ncp,Ncm,Mc,~]...
+ = PCMfilteringF(fCo,fCp,fCm,OTFo,OBJparaA,kC);
+
+%% doubling Fourier domain size if necessary
+OTFo = OTFdoubling(OTFo,DoubleMatSize);
+
+%% merging all 9 frequency components using generalized Wiener Filter
+[Fsum,Fperi,Fcent] = MergingHeptaletsF(fAof,fApf,fAmf,...
+ fBof,fBpf,fBmf,fCof,fCpf,fCmf,...
+ Ma,Mb,Mc,Nao,Nap,Nam,Nbo,Nbp,Nbm,...
+ Nco,Ncp,Ncm,kA,kB,kC,OBJparaA,OTFo);
+
+% Plotting SIM results
+SIMplot(Fsum,Fperi,Fcent,OTFo,kA,kB,kC,S1aTnoisy);
+break
+
+%% recontructed SIM images
+Dcent = real( ifft2(fftshift(Fcent)) );
+Dsum = real( ifft2(fftshift(Fsum)) );
+Dperi = real( ifft2(fftshift(Fperi)) );
+
+%{
+figure;
+mesh(log(abs(Fsum)))
+figure;
+mesh(log(abs(Fperi)))
+%}
+
+w = size(OTFo,1);
+h = 30; % pixel width to be trimed of the image edge to trim of artifacts
+figure;
+imshow(Dcent(h+1:t-h,h+1:t-h),[])
+title('Wiener Filtered wide-field')
+figure;
+imshow(Dsum(h+1:t-h,h+1:t-h),[])
+title('SIM image')
+figure;
+imshow(Dperi(h+1:t-h,h+1:t-h),[])
+title('SIM image (using only off-center frequency components)')
+
+%% appodizing the merged frequency components
+Index = 0.4;
+Kotf = OTFedgeF(OTFo);
+[FsumA] = ApodizationFunction(Fsum,kA,kB,kC,Kotf,Index);
+[FperiA] = ApodizationFunction(Fperi,kA,kB,kC,Kotf,Index);
+DsumA = real( ifft2(fftshift(FsumA)) );
+DperiA = real( ifft2(fftshift(FperiA)) );
+
+figure;
+imshow(DsumA(h+1:t-h,h+1:t-h),[])
+title('appodized SIM image')
+figure;
+imshow(DperiA(h+1:t-h,h+1:t-h),[])
+title('appodized SIM image (using only off-center frequency components)')
+break
+%% for writing the image files into *.tiff files
+%{
+h = 1*30;
+imwrite(im2uint16(mat2gray(Dcent(h+1:t-h,h+1:t-h))),['Dcent.tif'],'tiff');
+imwrite(im2uint16(mat2gray(Dsum(h+1:t-h,h+1:t-h))),['Dsum.tif'],'tiff');
+imwrite(im2uint16(mat2gray(Dperi(h+1:t-h,h+1:t-h))),['Dperi.tif'],'tiff');
+imwrite(im2uint16(mat2gray(DsumA(h+1:t-h,h+1:t-h))),['DsumA.tif'],'tiff');
+imwrite(im2uint16(mat2gray(DperiA(h+1:t-h,h+1:t-h))),['DperiA.tif'],'tiff');
+%}
+
+%% ploting FT of images
+fDIo = fftshift(fft2(DIo));
+fDIoT = fftshift(fft2(DIoT));
+fDIoTnoisy = fftshift(fft2(DIoTnoisy));
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+
+p = 10;
+minL1 = min(min( abs(fDIo).^(1/p) ));
+minL2 = min(min( abs(fDIoT).^(1/p) ));
+minL3 = min(min( abs(fDIoTnoisy).^(1/p) ));
+minL4 = min(min( abs(fS1aTnoisy).^(1/p) ));
+minL5 = min(min( abs(Fcent).^(1/p) ));
+minL6 = min(min( abs(Fsum).^(1/p) ));
+minL7 = min(min( abs(Fperi).^(1/p) ));
+maxL1 = max(max( abs(fDIo).^(1/p) ));
+maxL2 = max(max( abs(fDIoT).^(1/p) ));
+maxL3 = max(max( abs(fDIoTnoisy).^(1/p) ));
+maxL4 = max(max( abs(fS1aTnoisy).^(1/p) ));
+maxL5 = max(max( abs(Fcent).^(1/p) ));
+maxL6 = max(max( abs(Fsum).^(1/p) ));
+maxL7 = max(max( abs(Fperi).^(1/p) ));
+minL = min([minL1,minL2,minL3,minL4,minL5,minL6,minL7]);
+maxL = max([maxL1,maxL2,maxL3,maxL4,maxL5,maxL6,maxL7]);
+
+figure;
+imshow(abs(fDIo).^(1/p),[minL maxL])
+title('fDIo')
+figure;
+imshow(abs(fDIoT).^(1/p),[minL maxL])
+title('fDIoT')
+figure;
+imshow(abs(fDIoTnoisy).^(1/p),[minL maxL])
+title('fDIoTnoisy')
+figure;
+imshow(abs(fS1aTnoisy).^(1/p),[minL maxL])
+title('fS1aTnoisy')
+
+figure;
+imshow(abs(Fcent).^(1/p),[minL maxL])
+title('Fcent')
+figure;
+imshow(abs(Fsum).^(1/p),[minL maxL])
+title('Fsum')
+figure;
+imshow(abs(Fperi).^(1/p),[minL maxL])
+title('Fperi')
+figure;
+imshow(abs(FsumA).^(1/p),[minL maxL])
+title('FsumA')
+figure;
+imshow(abs(FperiA).^(1/p),[minL maxL])
+title('FperiA')
+
+%% for writing the image files into *.tiff files
+%{
+imwrite(im2uint16(mat2gray(abs(Fcent).^(1/p))),['Fcent.tif'],'tiff');
+imwrite(im2uint16(mat2gray(abs(Fsum).^(1/p))),['Fsum.tif'],'tiff');
+imwrite(im2uint16(mat2gray(abs(Fperi).^(1/p))),['Fperi.tif'],'tiff');
+imwrite(im2uint16(mat2gray(abs(FsumA).^(1/p))),['FsumA.tif'],'tiff');
+imwrite(im2uint16(mat2gray(abs(FperiA).^(1/p))),['FperiA.tif'],'tiff');
+%}
diff --git a/SIMbasic/SIMplot.m b/SIMbasic/SIMplot.m
new file mode 100644
index 0000000..b83fd74
--- /dev/null
+++ b/SIMbasic/SIMplot.m
@@ -0,0 +1,84 @@
+function [] = SIMplot(Fsum,Fperi,Fcent,OTFo,kA,kB,kC,S1aTnoisy)
+
+%% recontructed SIM images
+Dcent = real( ifft2(fftshift(Fcent)) );
+Dsum = real( ifft2(fftshift(Fsum)) );
+Dperi = real( ifft2(fftshift(Fperi)) );
+
+t = size(OTFo,1);
+h = 1*30; % for removing the image edges
+figure;
+imshow(Dcent(h+1:t-h,h+1:t-h),[])
+title('Wiener-Filtered wide-field')
+figure;
+imshow(Dsum(h+1:t-h,h+1:t-h),[])
+title('SIM image')
+figure;
+imshow(Dperi(h+1:t-h,h+1:t-h),[])
+title('SIM image (using only off-center frequency components)')
+
+% appodizing the merged frequency components
+Index = 0.4;
+Kotf = OTFedgeF(OTFo);
+[FsumA] = ApodizationFunction(Fsum,kA,kB,kC,Kotf,Index);
+[FperiA] = ApodizationFunction(Fperi,kA,kB,kC,Kotf,Index);
+DsumA = real( ifft2(fftshift(FsumA)) );
+DperiA = real( ifft2(fftshift(FperiA)) );
+
+figure;
+imshow(DsumA(h+1:t-h,h+1:t-h),[])
+title('appodized SIM image')
+figure;
+imshow(DperiA(h+1:t-h,h+1:t-h),[])
+title('appodized SIM image (using only off-center frequency components)')
+
+
+%% Frequency Plots
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+w = size(OTFo,1);
+wo = w/2;
+w1 = size(fS1aTnoisy,1);
+if ( w > w1 )
+ fS1aTnoisy0 = zeros(w,w);
+ fS1aTnoisy0(wo-w1/2+1:wo+w1/2,wo-w1/2+1:wo+w1/2) = fS1aTnoisy;
+ clear fS1aTnoisy
+ fS1aTnoisy = fS1aTnoisy0;
+ clear fS1aTnoisy0;
+end
+
+
+p = 10;
+minL1 = min(min( abs(fS1aTnoisy).^(1/p) ));
+minL2 = min(min( abs(Fcent).^(1/p) ));
+minL3 = min(min( abs(Fsum).^(1/p) ));
+minL4 = min(min( abs(Fperi).^(1/p) ));
+minL5 = min(min( abs(FsumA).^(1/p) ));
+maxL1 = max(max( abs(fS1aTnoisy).^(1/p) ));
+maxL2 = max(max( abs(Fcent).^(1/p) ));
+maxL3 = max(max( abs(Fsum).^(1/p) ));
+maxL4 = max(max( abs(Fperi).^(1/p) ));
+maxL5 = max(max( abs(FsumA).^(1/p) ));
+minL = min([minL1,minL2,minL3,minL4,minL5]);
+maxL = max([maxL1,maxL2,maxL3,maxL4,maxL5]);
+
+figure;
+imshow(abs(fS1aTnoisy).^(1/p),[minL maxL])
+title('fS1aTnoisy')
+
+figure;
+imshow(abs(Fcent).^(1/p),[minL maxL])
+title('Weiner Filtered frequency')
+figure;
+imshow(abs(Fsum).^(1/p),[minL maxL])
+title('SIM frequency')
+figure;
+imshow(abs(Fperi).^(1/p),[minL maxL])
+title('SIM (off-center frequency components)')
+figure;
+imshow(abs(FsumA).^(1/p),[minL maxL])
+title('appodized SIM frequency')
+figure;
+imshow(abs(FperiA).^(1/p),[minL maxL])
+title('appodized SIM (off-center frequency components)')
+
+
diff --git a/SIMbasic/SeparatedComponents2D.m b/SIMbasic/SeparatedComponents2D.m
new file mode 100644
index 0000000..14d3880
--- /dev/null
+++ b/SIMbasic/SeparatedComponents2D.m
@@ -0,0 +1,22 @@
+function [FiSMao,FiSMap,FiSMam] = SeparatedComponents2D(...
+ phaseShift,phaseShift0,FcS1aT,FcS2aT,FcS3aT)
+% Aim: Unmixing the frequency components of raw SIM images
+% phaseShift,phaseShift0: illumination phase shifts
+% FcS1aT,FcS2aT,FcS3aT: FT of raw SIM images
+% FiSMao,FiSMap,FiSMam: unmixed frequency components of raw SIM images
+
+phaseShift1 = phaseShift(1,1);
+phaseShift2 = phaseShift(2,1);
+MF = 1.0;
+%% Transformation Matrix
+M = 0.5*[1 0.5*MF*exp(-1i*phaseShift0) 0.5*MF*exp(+1i*phaseShift0);
+ 1 0.5*MF*exp(-1i*phaseShift1) 0.5*MF*exp(+1i*phaseShift1);
+ 1 0.5*MF*exp(-1i*phaseShift2) 0.5*MF*exp(+1i*phaseShift2)];
+
+%% Separting the components
+%===========================================================
+Minv = inv(M);
+
+FiSMao = Minv(1,1)*FcS1aT + Minv(1,2)*FcS2aT + Minv(1,3)*FcS3aT;
+FiSMap = Minv(2,1)*FcS1aT + Minv(2,2)*FcS2aT + Minv(2,3)*FcS3aT;
+FiSMam = Minv(3,1)*FcS1aT + Minv(3,2)*FcS2aT + Minv(3,3)*FcS3aT;
diff --git a/SIMbasic/TripletSNR0.m b/SIMbasic/TripletSNR0.m
new file mode 100644
index 0000000..6bf5b1d
--- /dev/null
+++ b/SIMbasic/TripletSNR0.m
@@ -0,0 +1,66 @@
+function [SIGao,SIGap2,SIGam2] = TripletSNR0(OBJpara,k2fa,OTFo,fDIp)
+
+% AIM: To obtain signal spectrums corresponding to central and off-center
+% frequency components
+% INPUT VARIABLES
+% OBJpara: object power parameters
+% k2fa: illumination frequency vector
+% OTFo: system OTF
+% fDIp: one of the off-center frequency component (utilized here only for
+% visual verification of computation)
+% OUTPUT VARIABLES
+% SIGao: signal spectrum corresponding to central frequency component
+% SIGap2,SIGam2: signal spectrums corresponding to off-center frequency components
+
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% object spectrum parameters
+Aobj = OBJpara(1);
+Bobj = OBJpara(2);
+
+% object spectrums (default)
+kv = k2fa(2) + 1i*k2fa(1); % vector along illumination direction
+Rp = abs(Cv-kv);
+Rm = abs(Cv+kv);
+OBJo = Aobj*(Ro.^Bobj);
+OBJp = Aobj*(Rp.^Bobj);
+OBJm = Aobj*(Rm.^Bobj);
+
+OBJo(wo+1,wo+1) = 0.25*OBJo(wo+2,wo+1) + 0.25*OBJo(wo+1,wo+2)...
+ + 0.25*OBJo(wo+0,wo+1) + 0.25*OBJo(wo+1,wo+0);
+
+k3 = round(k2fa);
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+OBJm(wo+1-k3(1),wo+1-k3(2)) = 0.25*OBJm(wo+2-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+2-k3(2))...
+ + 0.25*OBJm(wo+0-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+0-k3(2));
+
+% signal spectrums
+SIGao = OBJo.*OTFo;
+SIGap = OBJp.*OTFo;
+SIGam = OBJm.*OTFo;
+
+SIGap2 = circshift(SIGap,-k3);
+SIGam2 = circshift(SIGam,k3);
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot([0:w-1]-wo,log(abs(SIGap2(wo+pp,:))),'k--','LineWidth',3,'MarkerSize',6)
+plot([0:w-1]-wo,log(abs(fDIp(wo+pp,:))),'mo-','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/SIMbasic/WoFilterCenter.m b/SIMbasic/WoFilterCenter.m
new file mode 100644
index 0000000..8299473
--- /dev/null
+++ b/SIMbasic/WoFilterCenter.m
@@ -0,0 +1,57 @@
+function [FiSMaof,NoisePower] = WoFilterCenter(FiSMao,OTFo,co,OBJpara,SFo)
+% Aim: Wiener Filtering the central frequency component
+% INPUT VARIABLES
+% FiSMao: noisy central frequency component
+% OTFo: system OTF
+% co: Wiener filter constant [=1, for minimum RMS estimate]
+% SFo: scaling factor (not significant here, so set to 1)
+% OUTPUT VARIABLES
+% FiSMaof: Wiener Filtered estimate of FiSMao
+% NoisePower: avg. noise power in FiSMao
+
+w = size(FiSMao,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Ro = sqrt( (X-wo).^2 + (Y-wo).^2 );
+
+OTFpower = OTFo.*conj(OTFo);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = FiSMao.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Object Power determination
+Aobj = OBJpara(1);
+Bobj = OBJpara(2);
+Ro(wo+1,wo+1) = 1;
+OBJpower = Aobj*(Ro.^Bobj);
+OBJpower = OBJpower.^2;
+
+%% Wiener Filtering
+FiSMaof = FiSMao.*(SFo.*conj(OTFo)./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+
+%% for cross-checking filtered estimate visually
+%{
+WFilter = (SFo.^2.*OTFpower./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+FiSMao1 = FiSMaof.*OTFo;
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(FiSMaof(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao(wo+pp,:))),'o--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao1(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,10*WFilter(wo+pp,:),'k*--','LineWidth',2,'MarkerSize',4)
+title('WoFilter')
+legend('FiSMaof','FiSMao','FiSMao1')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/SIMbasic/WoFilterSideLobe.m b/SIMbasic/WoFilterSideLobe.m
new file mode 100644
index 0000000..fb717c5
--- /dev/null
+++ b/SIMbasic/WoFilterSideLobe.m
@@ -0,0 +1,54 @@
+function [FiSMaof,NoisePower] = WoFilterSideLobe(FiSMao,OTFo,co,OBJsideP,SFo)
+% Aim: Wiener Filtering the off-center frequency component
+% INPUT VARIABLES
+% FiSMao: noisy off-center frequency component
+% OTFo: system OTF
+% co: Wiener filter constant [=1, for minimum RMS estimate]
+% OBJsideP: avg. object spectrum
+% SFo: modulation factor
+% OUTPUT VARIABLES
+% FiSMaof: Wiener Filtered estimate of FiSMao
+% NoisePower: avg. noise power in FiSMao
+
+w = size(FiSMao,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Ro = sqrt( (X-wo).^2 + (Y-wo).^2 );
+
+OTFpower = OTFo.*conj(OTFo);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = FiSMao.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Object Power determination
+OBJpower = OBJsideP.^2;
+
+%% Wiener Filtering
+FiSMaof = FiSMao.*(SFo.*conj(OTFo)./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+
+%% for cross-checking filtered estimate visually
+%{
+WFilter = (SFo.^2.*OTFpower./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+FiSMao1 = FiSMaof.*OTFo;
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(FiSMaof(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao(wo+pp,:))),'o--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao1(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,10*WFilter(wo+pp,:),'k*--','LineWidth',2,'MarkerSize',4)
+title('WoFilter')
+legend('FiSMaof','FiSMao','FiSMao1')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/SIMbasic/testpat.tiff b/SIMbasic/testpat.tiff
new file mode 100644
index 0000000..84f1d62
Binary files /dev/null and b/SIMbasic/testpat.tiff differ
diff --git a/SIMexpt/ApodizationFunction.m b/SIMexpt/ApodizationFunction.m
new file mode 100644
index 0000000..172f322
--- /dev/null
+++ b/SIMexpt/ApodizationFunction.m
@@ -0,0 +1,30 @@
+function [FsumA] = ApodizationFunction(Fsum,k2fa,k2fb,k2fc,Kotf,Index)
+
+t = size(Fsum,1);
+
+OTFo_mask = OTFmaskShifted(0.*k2fa,Kotf,t);
+OTFap_mask = OTFmaskShifted(k2fa,Kotf,t);
+OTFam_mask = OTFmaskShifted(-k2fa,Kotf,t);
+OTFbp_mask = OTFmaskShifted(k2fb,Kotf,t);
+OTFbm_mask = OTFmaskShifted(-k2fb,Kotf,t);
+OTFcp_mask = OTFmaskShifted(k2fc,Kotf,t);
+OTFcm_mask = OTFmaskShifted(-k2fc,Kotf,t);
+ApoMask = OTFo_mask.*OTFap_mask.*OTFam_mask.*OTFbp_mask.*OTFbm_mask.*OTFcp_mask.*OTFcm_mask;
+clear OTFoz OTFo_mask OTFap_mask OTFam_mask OTFbp_mask OTFbm_mask OTFcp_mask OTFcm_mask
+%{
+figure;
+imshow(ApoMask,[ ])
+%}
+
+DistApoMask = bwdist(ApoMask);
+maxApoMask = max(max(DistApoMask));
+ApoFunc = double(DistApoMask./maxApoMask).^Index;
+%{
+figure;
+imshow(DistApoMask,[ ])
+figure;
+mesh(double(ApoFunc))
+%}
+
+FsumA = Fsum.*ApoFunc;
+clear ApoMask DistApoMask ApoFunc
\ No newline at end of file
diff --git a/SIMexpt/ApproxFreqDuplex.m b/SIMexpt/ApproxFreqDuplex.m
new file mode 100644
index 0000000..a912bd5
--- /dev/null
+++ b/SIMexpt/ApproxFreqDuplex.m
@@ -0,0 +1,28 @@
+function [maxK2,Ix,Iy] = ApproxFreqDuplex(FiSMap,Kotf)
+% AIM: approx. illumination frequency vector determination
+% INPUT VARIABLES
+% FiSMap: FT of raw SIM image
+% Kotf: OTF cut-off frequency
+% OUTPUT VARIABLES
+% maxK2: illumination frequency vector (approx)
+% Ix,Iy: coordinates of illumination frequency peaks
+
+FiSMap = abs(FiSMap);
+
+w = size(FiSMap,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+Ro = sqrt( (X-wo).^2 + (Y-wo).^2 );
+Z0 = Ro > round(0.5*Kotf);
+Z1 = X > wo;
+
+FiSMap = FiSMap.*Z0.*Z1;
+dumY = max( FiSMap,[],1 );
+[dummy Iy] = max(dumY);
+dumX = max( FiSMap,[],2 );
+[dummy Ix] = max(dumX);
+
+maxK2 = [Ix-(wo+1) Iy-(wo+1)];
diff --git a/SIMexpt/BgNormalization0.m b/SIMexpt/BgNormalization0.m
new file mode 100644
index 0000000..7d550b1
--- /dev/null
+++ b/SIMexpt/BgNormalization0.m
@@ -0,0 +1,39 @@
+function [S1aTnoisyA, S2aTnoisyA, S3aTnoisyA, ...
+ S1bTnoisyA, S2bTnoisyA, S3bTnoisyA, ...
+ S1cTnoisyA, S2cTnoisyA, S3cTnoisyA] = BgNormalization0(...
+ S1aTnoisy, S2aTnoisy, S3aTnoisy, ...
+ S1bTnoisy, S2bTnoisy, S3bTnoisy, ...
+ S1cTnoisy, S2cTnoisy, S3cTnoisy)
+
+m1a = mean2(S1aTnoisy);
+m2a = mean2(S2aTnoisy);
+m3a = mean2(S3aTnoisy);
+m1b = mean2(S1bTnoisy);
+m2b = mean2(S2bTnoisy);
+m3b = mean2(S3bTnoisy);
+m1c = mean2(S1cTnoisy);
+m2c = mean2(S2cTnoisy);
+m3c = mean2(S3cTnoisy);
+
+s1a = std(S1aTnoisy(:));
+s2a = std(S2aTnoisy(:));
+s3a = std(S3aTnoisy(:));
+s1b = std(S1bTnoisy(:));
+s2b = std(S2bTnoisy(:));
+s3b = std(S3bTnoisy(:));
+s1c = std(S1cTnoisy(:));
+s2c = std(S2cTnoisy(:));
+s3c = std(S3cTnoisy(:));
+
+mTemp = max([m1a m2a m3a m1b m2b m3b m1c m2c m3c]);
+sTemp = max([s1a s2a s3a s1b s2b s3b s1c s2c s3c]);
+
+S1aTnoisyA = (S1aTnoisy - m1a).*(sTemp./s1a) + mTemp;
+S2aTnoisyA = (S2aTnoisy - m2a).*(sTemp./s2a) + mTemp;
+S3aTnoisyA = (S3aTnoisy - m3a).*(sTemp./s3a) + mTemp;
+S1bTnoisyA = (S1bTnoisy - m1b).*(sTemp./s1b) + mTemp;
+S2bTnoisyA = (S2bTnoisy - m2b).*(sTemp./s2b) + mTemp;
+S3bTnoisyA = (S3bTnoisy - m3b).*(sTemp./s3b) + mTemp;
+S1cTnoisyA = (S1cTnoisy - m1c).*(sTemp./s1c) + mTemp;
+S2cTnoisyA = (S2cTnoisy - m2c).*(sTemp./s2c) + mTemp;
+S3cTnoisyA = (S3cTnoisy - m3c).*(sTemp./s3c) + mTemp;
\ No newline at end of file
diff --git a/SIMexpt/IlluminationFreqF.m b/SIMexpt/IlluminationFreqF.m
new file mode 100644
index 0000000..5b87d16
--- /dev/null
+++ b/SIMexpt/IlluminationFreqF.m
@@ -0,0 +1,38 @@
+function [k2fa] = IlluminationFreqF(S1aTnoisy,OTFo)
+% AIM: illumination frequency vector determination
+% INPUT VARIABLES
+% S1aTnoisy: raw SIM image
+% OTFo: system OTF
+% OUTPUT VARIABLE
+% k2fa: illumination frequency vector
+
+w = size(OTFo,1);
+wo = w/2;
+
+% computing PSFe for edge tapering SIM images (heuristically determined)
+PSFd = real(fftshift( ifft2(fftshift(OTFo.^10)) ));
+PSFd = PSFd/max(max(PSFd));
+PSFd = PSFd/sum(sum(PSFd));
+h = 30;
+PSFe = PSFd(wo-h+1:wo+h,wo-h+1:wo+h);
+
+% edge tapering raw SIM image
+S1aTnoisy_et = edgetaper(S1aTnoisy,PSFe);
+fS1aTnoisy_et = fftshift(fft2(S1aTnoisy_et));
+
+% OTF cut-off freq
+Kotf = OTFedgeF(OTFo);
+
+% Approx illumination frequency vector
+[k2fa,~,~] = ApproxFreqDuplex(fS1aTnoisy_et,Kotf);
+
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+% illumination frequency vector determination by optimizing
+% autocorrelation of fS1aTnoisy
+OPT = 1;
+PhaseKai2opt0 = @(k2fa0)PhaseKai2opt(k2fa0,fS1aTnoisy,OTFo,OPT);
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+k2fa0 = k2fa;
+[k2fa,fval] = fminsearch(PhaseKai2opt0,k2fa0,options);
+% k2a = sqrt(k2fa*k2fa')
diff --git a/SIMexpt/IlluminationPhaseF.m b/SIMexpt/IlluminationPhaseF.m
new file mode 100644
index 0000000..2898568
--- /dev/null
+++ b/SIMexpt/IlluminationPhaseF.m
@@ -0,0 +1,14 @@
+function [phaseA1] = IlluminationPhaseF(S1aTnoisy,k2fa)
+% AIM: illumination phase shift determination
+% INPUT VARIABLES
+% S1aTnoisy: raw SIM image
+% k2fa: illumination frequency vector
+% OUTPUT VARIABLE
+% phaseA1: illumination phase shift determined
+
+PatternPhaseOpt0 = @(phaseA0)PatternPhaseOpt(phaseA0,S1aTnoisy,k2fa);
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+phaseA0 = 0;
+[phaseA1,fval] = fminsearch(PatternPhaseOpt0,phaseA0,options);
+% phaseA1*180/pi
\ No newline at end of file
diff --git a/SIMexpt/MergingHeptaletsF.m b/SIMexpt/MergingHeptaletsF.m
new file mode 100644
index 0000000..d1b915b
--- /dev/null
+++ b/SIMexpt/MergingHeptaletsF.m
@@ -0,0 +1,65 @@
+function [Fsum,Fperi,Fcent] = MergingHeptaletsF(fDIo,fDIp,fDIm,fDBo,fDBp,fDBm,...
+ fDCo,fDCp,fDCm,Ma,Mb,Mc,npDIo,npDIp,npDIm,...
+ npDBo,npDBp,npDBm,npDCo,npDCp,npDCm,k2fa,k2fb,k2fc,OBJpara,OTFo)
+
+% AIM: To merge all 9 frequency components into one using generalized
+% Wiener Filter
+% INPUT VARIABLES
+% [fDIo fDIp fDIm
+% fDBo fDBp fDBm
+% fDCo fDCp fDCm]: nine frequency components
+% Ma,Mb,Mc: modulation factors for the three illumination orientations
+% [npDIo npDIp npDIm
+% npDBo npDBp npDBm
+% npDCo,npDCp,npDCm]: noise powers corresponding to nine frequency
+% components
+% k2fa,k2fb,k2fc: illumination frequency vectors for the three
+% illumination orientations
+% OBJpara: Object spectrum parameters
+% OTFo: system OTF
+% OUTPUT VARIABLES
+% Fsum: all nine frequency components merged into one using
+% generalised Wiener Filter
+% Fperi: six off-center frequency components merged into one using
+% generalised Wiener Filter
+% Fcent: averaged of the three central frequency components
+
+% obtain signal spectrums corresponding to central and off-center
+% frequency components
+[SIGao,SIGam2,SIGap2] = TripletSNR0(OBJpara,k2fa,OTFo,fDIm);
+[SIGbo,SIGbm2,SIGbp2] = TripletSNR0(OBJpara,k2fb,OTFo,fDBm);
+[SIGco,SIGcm2,SIGcp2] = TripletSNR0(OBJpara,k2fc,OTFo,fDCm);
+
+SIGap2 = Ma*SIGap2;
+SIGam2 = Ma*SIGam2;
+SIGbp2 = Mb*SIGbp2;
+SIGbm2 = Mb*SIGbm2;
+SIGcp2 = Mc*SIGcp2;
+SIGcm2 = Mc*SIGcm2;
+
+%% Generalized Wiener-Filter computation
+SNRao = SIGao.*conj(SIGao)./npDIo;
+SNRap = SIGap2.*conj(SIGap2)./npDIp;
+SNRam = SIGam2.*conj(SIGam2)./npDIm;
+
+SNRbo = SIGbo.*conj(SIGbo)./npDBo;
+SNRbp = SIGbp2.*conj(SIGbp2)./npDBp;
+SNRbm = SIGbm2.*conj(SIGbm2)./npDBm;
+
+SNRco = SIGco.*conj(SIGco)./npDCo;
+SNRcp = SIGcp2.*conj(SIGcp2)./npDCp;
+SNRcm = SIGcm2.*conj(SIGcm2)./npDCm;
+
+ComDeno = 0.01 + ( SNRao + SNRap + SNRam + SNRbo + SNRbp + SNRbm + SNRco + SNRcp + SNRcm );
+Fsum = fDIo.*SNRao + fDIp.*SNRap + fDIm.*SNRam...
+ + fDBo.*SNRbo + fDBp.*SNRbp + fDBm.*SNRbm...
+ + fDCo.*SNRco + fDCp.*SNRcp + fDCm.*SNRcm;
+Fsum = Fsum./ComDeno;
+
+ComPeri = 0.01 + ( SNRap + SNRam + SNRbp + SNRbm + SNRcp + SNRcm );
+Fperi = fDIp.*SNRap + fDIm.*SNRam + fDBp.*SNRbp + fDBm.*SNRbm...
+ + fDCp.*SNRcp + fDCm.*SNRcm;
+Fperi = Fperi./ComPeri;
+
+% averaged central frequency component
+Fcent = (fDIo+fDBo+fDCo)/3;
diff --git a/SIMexpt/ModulationFactor.m b/SIMexpt/ModulationFactor.m
new file mode 100644
index 0000000..4ce0313
--- /dev/null
+++ b/SIMexpt/ModulationFactor.m
@@ -0,0 +1,86 @@
+function Mm = ModulationFactor(fDp,k1a,OBJparaA,OTFo)
+% AIM: Determination of modulation factor
+% INPUT VARIABLES
+% fDp: off-center frequency component
+% k1a: illumination frequency vector
+% OBJparaA: Object power parameters
+% OTFo: system OTF
+% OUTPUT VARIABLE
+% Mm: modulation factor
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% magnitude of illumination vector
+k2 = sqrt(k1a*k1a');
+
+% vector along illumination direction
+kv = k1a(2) + 1i*k1a(1);
+Rp = abs(Cv+kv);
+
+% Object power parameters
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+
+% Object spectrum
+OBJp = Aobj*(Rp+0).^Bobj;
+
+% illumination vector rounded to nearest pixel
+k3 = -round(k1a);
+
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+
+% signal spectrum
+SIGap = OBJp.*OTFo;
+
+% figure;
+% mesh(log(abs(OBJp)))
+% figure;
+% mesh(log(abs(fDp)))
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = fDp.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Noise free object power computation
+Fpower = fDp.*conj(fDp) - NoisePower;
+fDp = sqrt(abs(Fpower));
+
+% frequency range over which signal power matching is done to estimate
+% modulation factor
+Zmask = (Ro > 0.2*k2).*(Ro < 0.8*k2).*(Rp > 0.2*k2);
+
+% least square approximation for modulation factor
+Mm = sum(sum(SIGap.*abs(fDp).*Zmask));
+Mm = Mm./sum(sum(SIGap.^2.*Zmask));
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(fDp(wo+pp,:))),'o-','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(SIGap(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,log(abs(Mm*SIGap(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+grid on
+box on
+kk
+%}
+
+
+
diff --git a/SIMexpt/OBJparaOpt.m b/SIMexpt/OBJparaOpt.m
new file mode 100644
index 0000000..8c17130
--- /dev/null
+++ b/SIMexpt/OBJparaOpt.m
@@ -0,0 +1,48 @@
+function Esum = OBJparaOpt(OBJpara0,fDIoTnoisy,OTFo)
+% AIM: Determined Sum of Squared Errors (SSE) between `actual signal power' and
+% `approximated signal power'
+% INPUT VARIABLES
+% OBJpara0: [Aobj Bobj], Object power parameters
+% fDIoTnoisy: FT of central frequency component
+% OTFo: system OTF
+% OUTPUT VARIABLE
+% Esum: SSE between `actual signal spectrum' and `approximated signal
+% spectrum'
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+Ro(wo+1,wo+1) = 1; % to avoid nan
+
+% approximated signal power calculation
+Aobj = OBJpara0(1);
+Bobj = OBJpara0(2);
+OBJpower = Aobj*(Ro.^Bobj);
+SIGpower = OBJpower.*OTFo;
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% range of frequency over which SSE is computed
+Zloop = (Ro<0.75*Kotf).*(Ro>0.25*Kotf);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = fDIoTnoisy.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Noise free object power computation
+Fpower = fDIoTnoisy.*conj(fDIoTnoisy) - NoisePower;
+fDIoTnoisy = sqrt(abs(Fpower));
+
+% SSE computation
+Error = fDIoTnoisy - SIGpower;
+Esum = sum(sum((Error.^2./Ro).*Zloop));
+
diff --git a/SIMexpt/OBJpowerPara.m b/SIMexpt/OBJpowerPara.m
new file mode 100644
index 0000000..9ef86ce
--- /dev/null
+++ b/SIMexpt/OBJpowerPara.m
@@ -0,0 +1,48 @@
+function OBJparaA = OBJpowerPara(fDIoTnoisy,OTFo)
+% AIM: determination of object power parameters Aobj and Bobj
+% INPUT VARIABLES
+% fDIoTnoisy: FT of central frequency component
+% OTFo: system OTFo
+% OUTPUT VARIABLES
+% OBJparaA: [Aobj Bobj]
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+%% object power parameters through optimization
+OBJparaOpt0 = @(OBJpara0)OBJparaOpt(OBJpara0,fDIoTnoisy,OTFo);
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+
+% obtaining crude initial guesses for Aobj and Bobj
+Zm = (Ro>0.3*Kotf).*(Ro<0.4*Kotf);
+Aobj = sum(sum(abs(fDIoTnoisy.*Zm)))./sum(sum(Zm));
+Bobj = -0.5;
+OBJpara0 = [Aobj Bobj];
+
+% optimization step
+[OBJparaA,fval] = fminsearch(OBJparaOpt0,OBJpara0,options);
+
+%% plot for cross-checking the result
+%{
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+OBJo = Aobj*(Ro.^Bobj);
+SIGp = OBJo.*OTFo;
+pp = 3;
+figure;
+hold on
+plot([0:w-1]-wo,log(abs(fDIoTnoisy(wo+pp,:))),'k--','LineWidth',3,'MarkerSize',6)
+plot([0:w-1]-wo,log(abs(SIGp(wo+pp,:))),'mo-','LineWidth',2,'MarkerSize',6)
+legend('acutal signal power','avg. signal power')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/SIMexpt/OTF.tif b/SIMexpt/OTF.tif
new file mode 100644
index 0000000..424998d
Binary files /dev/null and b/SIMexpt/OTF.tif differ
diff --git a/SIMexpt/OTFdoubling.m b/SIMexpt/OTFdoubling.m
new file mode 100644
index 0000000..58aaa6e
--- /dev/null
+++ b/SIMexpt/OTFdoubling.m
@@ -0,0 +1,17 @@
+function OTF2 = OTFdoubling(OTFo,DoubleMatSize)
+
+%% embeds OTFo in doubled range of frequencies
+
+w = size(OTFo,1);
+wo = w/2;
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ OTF2 = zeros(2*w,2*w);
+ OTF2(wo+1:w+wo,wo+1:w+wo) = OTFo;
+ clear OTFo
+else
+ OTF2 = OTFo;
+ clear OTFo
+end
\ No newline at end of file
diff --git a/SIMexpt/OTFedgeF.m b/SIMexpt/OTFedgeF.m
new file mode 100644
index 0000000..8148030
--- /dev/null
+++ b/SIMexpt/OTFedgeF.m
@@ -0,0 +1,13 @@
+function Kotf = OTFedgeF(OTFo)
+
+w = size(OTFo,1);
+wo = w/2;
+
+OTF1 = OTFo(wo+1,:);
+OTFmax = max(max(abs(OTFo)));
+OTFtruncate = 0.02;
+i = 1;
+while ( abs(OTF1(1,i))<OTFtruncate*OTFmax )
+ Kotf = wo+1-i;
+ i = i + 1;
+end
\ No newline at end of file
diff --git a/SIMexpt/OTFmaskShifted.m b/SIMexpt/OTFmaskShifted.m
new file mode 100644
index 0000000..a60ac98
--- /dev/null
+++ b/SIMexpt/OTFmaskShifted.m
@@ -0,0 +1,13 @@
+function [OTFap_mask] = OTFmaskShifted(k2fa,Kotf,w)
+% OTFap_mask = complement of shifted OTF
+
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+Iy = k2fa(1,1) + (wo+1);
+Ix = k2fa(1,2) + (wo+1);
+
+Ro = sqrt( (X+1-Ix).^2 + (Y+1-Iy).^2 );
+OTFap_mask = Ro > Kotf;
\ No newline at end of file
diff --git a/SIMexpt/OTFpost.m b/SIMexpt/OTFpost.m
new file mode 100644
index 0000000..2498868
--- /dev/null
+++ b/SIMexpt/OTFpost.m
@@ -0,0 +1,20 @@
+function OTFo = OTFpost(OTFo)
+%% AIM: To flaten out the peripheral rippling frequencies
+
+OTFo = OTFo./max(OTFo(:));
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+Zmask = (Ro>(Kotf+20)).*(Ro<wo);
+OTFperi = sum(sum(OTFo.*Zmask))./sum(sum(Zmask));
+Zperi = (Ro>(wo-10));
+OTFo = OTFo.*(1-Zperi) + OTFperi.*Zperi;
\ No newline at end of file
diff --git a/SIMexpt/PCMfilteringF.m b/SIMexpt/PCMfilteringF.m
new file mode 100644
index 0000000..2e3fbd4
--- /dev/null
+++ b/SIMexpt/PCMfilteringF.m
@@ -0,0 +1,157 @@
+function [fDof,fDp2,fDm2,npDo,npDp,npDm,Mm,DoubleMatSize]...
+ = PCMfilteringF(fDo,fDp,fDm,OTFo,OBJparaA,kA)
+
+% AIM: obtaining Wiener Filtered estimates of noisy frequency components
+% INPUT VARIABLES
+% fDo,fDp,fDm: noisy estimates of separated frequency components
+% OTFo: system OTF
+% OBJparaA: object power parameters
+% kA: illumination vector
+% OUTPUT VARIABLES
+% fDof,fDp2,fDm2: Wiener Filtered estimates of fDo,fDp,fDm
+% npDo,npDp,npDm: avg. noise power in fDo,fDp,fDm
+% Mm: illumination modulation factor
+% DoubleMatSize: parameter for doubling FT size if necessary
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% suppressing out-of-focus signal of off-center components using
+% G is a notch-filter (determined heuristically)
+G = 1 - exp(-0.05*Ro.^1.2);
+fDp = fDp.*G;
+fDm = fDm.*G;
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+
+% Wiener Filtering central frequency component
+SFo = 1;
+co = 1.0;
+[fDof,npDo] = WoFilterCenter(fDo,OTFo,co,OBJparaA,SFo);
+
+% modulation factor determination
+Mm = ModulationFactor(fDp,kA,OBJparaA,OTFo);
+
+%% Duplex power (default)
+kv = kA(2) + 1i*kA(1); % vector along illumination direction
+Rp = abs(Cv-kv);
+Rm = abs(Cv+kv);
+OBJp = Aobj*(Rp.^Bobj);
+OBJm = Aobj*(Rm.^Bobj);
+k3 = round(kA);
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+OBJm(wo+1-k3(1),wo+1-k3(2)) = 0.25*OBJm(wo+2-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+2-k3(2))...
+ + 0.25*OBJm(wo+0-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+0-k3(2));
+
+% Filtering side lobes (off-center frequency components)
+SFo = Mm;
+[fDpf,npDp] = WoFilterSideLobe(fDp,OTFo,co,OBJm,SFo);
+[fDmf,npDm] = WoFilterSideLobe(fDm,OTFo,co,OBJp,SFo);
+
+%% doubling Fourier domain size if necessary
+DoubleMatSize = 0;
+if ( 2*Kotf > wo )
+ DoubleMatSize = 1; % 1 for doubling fourier domain size, 0 for keeping it unchanged
+end
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ to = t/2;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ [U,V] = meshgrid(u,v);
+ fDoTemp = zeros(2*w,2*w);
+ fDpTemp = zeros(2*w,2*w);
+ fDmTemp = zeros(2*w,2*w);
+ OTFtemp = zeros(2*w,2*w);
+ fDoTemp(wo+1:w+wo,wo+1:w+wo) = fDof;
+ fDpTemp(wo+1:w+wo,wo+1:w+wo) = fDpf;
+ fDmTemp(wo+1:w+wo,wo+1:w+wo) = fDmf;
+ OTFtemp(wo+1:w+wo,wo+1:w+wo) = OTFo;
+ clear fDof fDpf fDmf OTFo w wo x y X Y
+ fDof = fDoTemp;
+ fDpf = fDpTemp;
+ fDmf = fDmTemp;
+ OTFo = OTFtemp;
+ clear fDoTemp fDpTemp fDmTemp OTFtemp
+else
+ t = w;
+ to = t/2;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ [U,V] = meshgrid(u,v);
+end
+
+% Shifting the off-center frequency components to their correct location
+fDp1 = fft2(ifft2(fDpf).*exp( +1i.*2*pi*(kA(2)/t.*(U-to) + kA(1)/t.*(V-to)) ));
+fDm1 = fft2(ifft2(fDmf).*exp( -1i.*2*pi*(kA(2)/t.*(U-to) + kA(1)/t.*(V-to)) ));
+%{
+figure;
+mesh(log(abs(fDp1)))
+
+pp = 3;
+figure;
+hold on
+plot(u-to,log(abs(fDof(to+pp,:))),'go-','LineWidth',2,'MarkerSize',6)
+plot(u-to,log(abs(fDp1(to+pp,:))),'o--','LineWidth',2,'MarkerSize',4)
+grid on
+box on
+% kk
+%}
+
+%% phase matching
+Cv = (U-to) + 1i*(V-to);
+Ro = abs(Cv);
+Rp = abs(Cv-kv);
+k2 = sqrt(kA*kA');
+
+% frequency range over which corrective phase is determined
+Zmask = (Ro < 0.8*k2).*(Rp < 0.8*k2);
+
+% corrective phase
+Angle0 = angle( sum(sum( fDof.*conj(fDp1).*Zmask )) );
+
+% phase correction
+fDp2 = exp(+1i*Angle0).*fDp1;
+fDm2 = exp(-1i*Angle0).*fDm1;
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot(u-to,angle(fDof(to+pp,:)).*180/pi,'o-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(fDp1(to+pp,:)).*180/pi,'ro-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(fDp2(to+pp,:)).*180/pi,'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+figure;
+hold on
+plot3(u-to,real(fDof(to+1,:)),imag(fDof(to+1,:)),'o-','LineWidth',2,'MarkerSize',8)
+plot3(u-to,real(fDp1(to+1,:)),imag(fDp1(to+1,:)),'ro-','LineWidth',2,'MarkerSize',8)
+plot3(u-to,real(fDp2(to+1,:)),imag(fDp2(to+1,:)),'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+
+ff1 = fDof.*conj(fDp1);
+ff2 = fDof.*conj(fDp2);
+figure;
+hold on
+plot(u-to,angle(ff1(to+pp,:)).*180/pi,'ro-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(ff2(to+pp,:)).*180/pi,'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/SIMexpt/PCMseparateF.m b/SIMexpt/PCMseparateF.m
new file mode 100644
index 0000000..2ec4c3d
--- /dev/null
+++ b/SIMexpt/PCMseparateF.m
@@ -0,0 +1,62 @@
+function [fDo,fDp,fDm,kA]...
+ = PCMseparateF(S1aTnoisy,S2aTnoisy,S3aTnoisy,OTFo)
+% AIM: obtaining the noisy estimates of three frequency components
+% INPUT VARIABLES
+% S1aTnoisy,S2aTnoisy,S3aTnoisy: 3 raw SIM images with identical
+% illumination pattern orientation
+% but different phase shifts
+% OTFo: system OTF
+% OUTPUT VARIABLES
+% fDo,fDp,fDm: noisy estimates of separated frequency components
+% kA: (averaged) illumination frequency vector
+
+w = size(S1aTnoisy,1);
+wo = w/2;
+
+%% Determination of illumination frequency vectors
+[k1a] = IlluminationFreqF(S1aTnoisy,OTFo);
+[k2a] = IlluminationFreqF(S2aTnoisy,OTFo);
+[k3a] = IlluminationFreqF(S3aTnoisy,OTFo);
+
+% mean illumination frequency vector
+kA = (k1a + k2a + k3a)/3;
+
+%% determination of illumination phase shifts
+[phase1A] = IlluminationPhaseF(S1aTnoisy,kA);
+[phase2A] = IlluminationPhaseF(S2aTnoisy,kA);
+[phase3A] = IlluminationPhaseF(S3aTnoisy,kA);
+% for display in command window
+phaseA = [phase1A; phase2A; phase3A];
+phaseA*180/pi
+
+% computing PSFe for edge tapering SIM images
+PSFd = real(fftshift( ifft2(fftshift(OTFo.^3)) ));
+PSFd = PSFd/max(max(PSFd));
+PSFd = PSFd/sum(sum(PSFd));
+h = 30;
+PSFe = PSFd(wo-h+1:wo+h,wo-h+1:wo+h);
+
+% edge tapering raw SIM images
+S1aTnoisy = edgetaper(S1aTnoisy,PSFe);
+S2aTnoisy = edgetaper(S2aTnoisy,PSFe);
+S3aTnoisy = edgetaper(S3aTnoisy,PSFe);
+
+
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+fS2aTnoisy = fftshift(fft2(S2aTnoisy));
+fS3aTnoisy = fftshift(fft2(S3aTnoisy));
+
+%% Separating the three frequency components
+phaseShift0 = phase1A;
+phaseShift = [phase2A; phase3A];
+[fDo,fDp,fDm] = SeparatedComponents2D(...
+ phaseShift,phaseShift0,fS1aTnoisy,fS2aTnoisy,fS3aTnoisy);
+%% for visual verification
+%{
+figure;
+mesh(log(abs(fDo)))
+figure;
+mesh(log(abs(fDp)))
+figure;
+mesh(log(abs(fDm)))
+%}
diff --git a/SIMexpt/PatternPhaseOpt.m b/SIMexpt/PatternPhaseOpt.m
new file mode 100644
index 0000000..0376834
--- /dev/null
+++ b/SIMexpt/PatternPhaseOpt.m
@@ -0,0 +1,13 @@
+function CCop = PatternPhaseOpt(phaseA,S1aTnoisy,k2fa)
+
+w = size(S1aTnoisy,1);
+wo = w/2;
+
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+sAo = cos( 2*pi*(k2fa(2).*(X-wo)+k2fa(1).*(Y-wo))./w + phaseA );
+% S1aTnoisy = real( ifft2(fftshift(fS1aTnoisy)) );
+S1aTnoisy = S1aTnoisy - mean2(S1aTnoisy);
+CCop = -sum(sum(S1aTnoisy.*sAo));
\ No newline at end of file
diff --git a/SIMexpt/PhaseKai2opt.m b/SIMexpt/PhaseKai2opt.m
new file mode 100644
index 0000000..1046bdf
--- /dev/null
+++ b/SIMexpt/PhaseKai2opt.m
@@ -0,0 +1,46 @@
+function [CCop] = PhaseKai2opt(k2fa,fS1aTnoisy,OTFo,OPT)
+% Aim: Compute autocorrelation of FT of raw SIM images
+% k2fa: illumination frequency vector
+% fS1aTnoisy: FT of raw SIM image
+% OTFo: system OTF
+% OPT: acronym for `OPTIMIZE'; to be set to 1 when this function is used
+% for optimization, or else to 0
+% CCop: autocorrelation of fS1aTnoisy
+
+w = size(fS1aTnoisy,1);
+wo = w/2;
+
+fS1aTnoisy = fS1aTnoisy.*(1-1*OTFo.^10);
+fS1aT = fS1aTnoisy.*conj(OTFo);
+
+Kotf = OTFedgeF(OTFo);
+DoubleMatSize = 0;
+if ( 2*Kotf > wo )
+ DoubleMatSize = 1; % 1 for doubling fourier domain size, 0 for keeping it unchanged
+end
+
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ fS1aT_temp = zeros(t,t);
+ fS1aT_temp(wo+1:w+wo,wo+1:w+wo) = fS1aT;
+ clear fS1aT
+ fS1aT = fS1aT_temp;
+ clear fS1aT_temp
+else
+ t = w;
+end
+to = t/2;
+u = linspace(0,t-1,t);
+v = linspace(0,t-1,t);
+[U,V] = meshgrid(u,v);
+S1aT = exp( -1i.*2*pi*( k2fa(2)/t.*(U-to)+k2fa(1)/t.*(V-to) ) ).*ifft2(fS1aT);
+fS1aT0 = fft2( S1aT );
+
+mA = sum(sum( fS1aT.*conj(fS1aT0) ));
+mA = mA./sum(sum( fS1aT0.*conj(fS1aT0) ));
+
+if (OPT > 0)
+ CCop = -abs(mA);
+else
+ CCop = mA;
+end
\ No newline at end of file
diff --git a/SIMexpt/SIMo.m b/SIMexpt/SIMo.m
new file mode 100644
index 0000000..a2aafd2
--- /dev/null
+++ b/SIMexpt/SIMo.m
@@ -0,0 +1,88 @@
+clear all
+close all
+clc
+
+%% Loading system OTF file
+OTFo = double(imread('OTF.tif'));
+OTFo = OTFpost(OTFo);
+
+%% Read Expt. Raw SIM Images
+aa1 = imread('sim01z4.tif');
+aa = aa1(:,:,2);
+bb = uint16(zeros(512,512,9));
+for ii=1:3
+ for jj=1:3
+ bb(1:512,1:512,(ii-1)*3+jj)=aa((ii-1)*512+1:ii*512,(jj-1)*512+1:jj*512,1);
+ end
+end
+
+% separating the raw SIM images
+S1aTnoisy = double( bb(:,:,1) );
+S2aTnoisy = double( bb(:,:,2) );
+S3aTnoisy = double( bb(:,:,3) );
+S1bTnoisy = double( bb(:,:,4) );
+S2bTnoisy = double( bb(:,:,5) );
+S3bTnoisy = double( bb(:,:,6) );
+S1cTnoisy = double( bb(:,:,7) );
+S2cTnoisy = double( bb(:,:,8) );
+S3cTnoisy = double( bb(:,:,9) );
+clear aa aa1 bb
+
+%% Intensity Normalization
+[S1aTnoisy, S2aTnoisy, S3aTnoisy, ...
+ S1bTnoisy, S2bTnoisy, S3bTnoisy, ...
+ S1cTnoisy, S2cTnoisy, S3cTnoisy] = BgNormalization0(...
+ S1aTnoisy, S2aTnoisy, S3aTnoisy, ...
+ S1bTnoisy, S2bTnoisy, S3bTnoisy, ...
+ S1cTnoisy, S2cTnoisy, S3cTnoisy);
+
+%% Background subtraction
+SE = strel('disk',20);
+S1aTnoisy = S1aTnoisy - imopen(S1aTnoisy,SE);
+S2aTnoisy = S2aTnoisy - imopen(S2aTnoisy,SE);
+S3aTnoisy = S3aTnoisy - imopen(S3aTnoisy,SE);
+S1bTnoisy = S1bTnoisy - imopen(S1bTnoisy,SE);
+S2bTnoisy = S2bTnoisy - imopen(S2bTnoisy,SE);
+S3bTnoisy = S3bTnoisy - imopen(S3bTnoisy,SE);
+S1cTnoisy = S1cTnoisy - imopen(S1cTnoisy,SE);
+S2cTnoisy = S2cTnoisy - imopen(S2cTnoisy,SE);
+S3cTnoisy = S3cTnoisy - imopen(S3cTnoisy,SE);
+
+
+%% obtaining the noisy estimates of three frequency components
+[fAo,fAp,fAm,kA]...
+ = PCMseparateF(S1aTnoisy,S2aTnoisy,S3aTnoisy,OTFo);
+[fBo,fBp,fBm,kB]...
+ = PCMseparateF(S1bTnoisy,S2bTnoisy,S3bTnoisy,OTFo);
+[fCo,fCp,fCm,kC]...
+ = PCMseparateF(S1cTnoisy,S2cTnoisy,S3cTnoisy,OTFo);
+
+% averaging the central frequency components
+fCent = (fAo + fBo + fCo)/3;
+
+% Object power parameters determination
+OBJparaA = OBJpowerPara(fCent,OTFo);
+
+%% Wiener Filtering the noisy frequency components
+[fAof,fApf,fAmf,Nao,Nap,Nam,Ma,DoubleMatSize]...
+ = PCMfilteringF(fAo,fAp,fAm,OTFo,OBJparaA,kA);
+[fBof,fBpf,fBmf,Nbo,Nbp,Nbm,Mb,~]...
+ = PCMfilteringF(fBo,fBp,fBm,OTFo,OBJparaA,kB);
+[fCof,fCpf,fCmf,Nco,Ncp,Ncm,Mc,~]...
+ = PCMfilteringF(fCo,fCp,fCm,OTFo,OBJparaA,kC);
+
+
+%% doubling Fourier domain size if necessary
+OTFo = OTFdoubling(OTFo,DoubleMatSize);
+
+%% merging all 9 frequency components using generalized Wiener Filter
+[Fsum,Fperi,Fcent] = MergingHeptaletsF(fAof,fApf,fAmf,...
+ fBof,fBpf,fBmf,fCof,fCpf,fCmf,...
+ Ma,Mb,Mc,Nao,Nap,Nam,Nbo,Nbp,Nbm,...
+ Nco,Ncp,Ncm,kA,kB,kC,OBJparaA,OTFo);
+
+
+% Plotting SIM results
+SIMplot(Fsum,Fperi,Fcent,OTFo,kA,kB,kC,S1aTnoisy);
+
+
diff --git a/SIMexpt/SIMplot.m b/SIMexpt/SIMplot.m
new file mode 100644
index 0000000..9a5f42f
--- /dev/null
+++ b/SIMexpt/SIMplot.m
@@ -0,0 +1,84 @@
+function [] = SIMplot(Fsum,Fperi,Fcent,OTFo,kA,kB,kC,S1aTnoisy)
+
+%% recontructed SIM images
+Dcent = real( ifft2(fftshift(Fcent)) );
+Dsum = real( ifft2(fftshift(Fsum)) );
+Dperi = real( ifft2(fftshift(Fperi)) );
+
+t = size(OTFo,1);
+h = 1*30;
+figure;
+imshow(Dcent(h+1:t-h,h+1:t-h),[])
+title('Wiener-Filtered wide-field')
+figure;
+imshow(Dsum(h+1:t-h,h+1:t-h),[])
+title('SIM image')
+figure;
+imshow(Dperi(h+1:t-h,h+1:t-h),[])
+title('SIM image (using only off-center frequency components)')
+
+% appodizing the merged frequency components
+Index = 0.4;
+Kotf = OTFedgeF(OTFo);
+[FsumA] = ApodizationFunction(Fsum,kA,kB,kC,Kotf,Index);
+[FperiA] = ApodizationFunction(Fperi,kA,kB,kC,Kotf,Index);
+DsumA = real( ifft2(fftshift(FsumA)) );
+DperiA = real( ifft2(fftshift(FperiA)) );
+
+figure;
+imshow(DsumA(h+1:t-h,h+1:t-h),[])
+title('appodized SIM image')
+figure;
+imshow(DperiA(h+1:t-h,h+1:t-h),[])
+title('appodized SIM image (using only off-center frequency components)')
+
+
+%% Frequency Plots
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+w = size(OTFo,1);
+wo = w/2;
+w1 = size(fS1aTnoisy,1);
+if ( w > w1 )
+ fS1aTnoisy0 = zeros(w,w);
+ fS1aTnoisy0(wo-w1/2+1:wo+w1/2,wo-w1/2+1:wo+w1/2) = fS1aTnoisy;
+ clear fS1aTnoisy
+ fS1aTnoisy = fS1aTnoisy0;
+ clear fS1aTnoisy0;
+end
+
+
+p = 10;
+minL1 = min(min( abs(fS1aTnoisy).^(1/p) ));
+minL2 = min(min( abs(Fcent).^(1/p) ));
+minL3 = min(min( abs(Fsum).^(1/p) ));
+minL4 = min(min( abs(Fperi).^(1/p) ));
+minL5 = min(min( abs(FsumA).^(1/p) ));
+maxL1 = max(max( abs(fS1aTnoisy).^(1/p) ));
+maxL2 = max(max( abs(Fcent).^(1/p) ));
+maxL3 = max(max( abs(Fsum).^(1/p) ));
+maxL4 = max(max( abs(Fperi).^(1/p) ));
+maxL5 = max(max( abs(FsumA).^(1/p) ));
+minL = min([minL1,minL2,minL3,minL4,minL5]);
+maxL = max([maxL1,maxL2,maxL3,maxL4,maxL5]);
+
+figure;
+imshow(abs(fS1aTnoisy).^(1/p),[minL maxL])
+title('fS1aTnoisy')
+
+figure;
+imshow(abs(Fcent).^(1/p),[minL maxL])
+title('Weiner Filtered frequency')
+figure;
+imshow(abs(Fsum).^(1/p),[minL maxL])
+title('SIM frequency')
+figure;
+imshow(abs(Fperi).^(1/p),[minL maxL])
+title('SIM (off-center frequency components)')
+figure;
+imshow(abs(FsumA).^(1/p),[minL maxL])
+title('appodized SIM frequency')
+figure;
+imshow(abs(FperiA).^(1/p),[minL maxL])
+title('appodized SIM (off-center frequency components)')
+
+
diff --git a/SIMexpt/SeparatedComponents2D.m b/SIMexpt/SeparatedComponents2D.m
new file mode 100644
index 0000000..14d3880
--- /dev/null
+++ b/SIMexpt/SeparatedComponents2D.m
@@ -0,0 +1,22 @@
+function [FiSMao,FiSMap,FiSMam] = SeparatedComponents2D(...
+ phaseShift,phaseShift0,FcS1aT,FcS2aT,FcS3aT)
+% Aim: Unmixing the frequency components of raw SIM images
+% phaseShift,phaseShift0: illumination phase shifts
+% FcS1aT,FcS2aT,FcS3aT: FT of raw SIM images
+% FiSMao,FiSMap,FiSMam: unmixed frequency components of raw SIM images
+
+phaseShift1 = phaseShift(1,1);
+phaseShift2 = phaseShift(2,1);
+MF = 1.0;
+%% Transformation Matrix
+M = 0.5*[1 0.5*MF*exp(-1i*phaseShift0) 0.5*MF*exp(+1i*phaseShift0);
+ 1 0.5*MF*exp(-1i*phaseShift1) 0.5*MF*exp(+1i*phaseShift1);
+ 1 0.5*MF*exp(-1i*phaseShift2) 0.5*MF*exp(+1i*phaseShift2)];
+
+%% Separting the components
+%===========================================================
+Minv = inv(M);
+
+FiSMao = Minv(1,1)*FcS1aT + Minv(1,2)*FcS2aT + Minv(1,3)*FcS3aT;
+FiSMap = Minv(2,1)*FcS1aT + Minv(2,2)*FcS2aT + Minv(2,3)*FcS3aT;
+FiSMam = Minv(3,1)*FcS1aT + Minv(3,2)*FcS2aT + Minv(3,3)*FcS3aT;
diff --git a/SIMexpt/TripletSNR0.m b/SIMexpt/TripletSNR0.m
new file mode 100644
index 0000000..5b6069e
--- /dev/null
+++ b/SIMexpt/TripletSNR0.m
@@ -0,0 +1,66 @@
+function [SIGao,SIGap2,SIGam2] = TripletSNR0(OBJpara,k2fa,OTFo,fDIp)
+
+% AIM: To obtain signal spectrums corresponding to central and off-center
+% frequency components
+% INPUT VARIABLES
+% OBJpara: object power parameters
+% k2fa: illumination frequency vector
+% OTFo: system OTF
+% fDIp: one of the off-center frequency component (utilized here only for
+% visual verification of computation)
+% OUTPUT VARIABLES
+% SIGao: signal spectrum corresponding to central frequency component
+% SIGap2,SIGam2: signal spectrums corresponding to off-center frequency
+% components
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% object spectrum parameters
+Aobj = OBJpara(1);
+Bobj = OBJpara(2);
+
+% object spectrums (default)
+kv = k2fa(2) + 1i*k2fa(1); % vector along illumination direction
+Rp = abs(Cv-kv);
+Rm = abs(Cv+kv);
+OBJo = Aobj*(Ro.^Bobj);
+OBJp = Aobj*(Rp.^Bobj);
+OBJm = Aobj*(Rm.^Bobj);
+
+OBJo(wo+1,wo+1) = 0.25*OBJo(wo+2,wo+1) + 0.25*OBJo(wo+1,wo+2)...
+ + 0.25*OBJo(wo+0,wo+1) + 0.25*OBJo(wo+1,wo+0);
+
+k3 = round(k2fa);
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+OBJm(wo+1-k3(1),wo+1-k3(2)) = 0.25*OBJm(wo+2-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+2-k3(2))...
+ + 0.25*OBJm(wo+0-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+0-k3(2));
+
+% signal spectrums
+SIGao = OBJo.*OTFo;
+SIGap = OBJp.*OTFo;
+SIGam = OBJm.*OTFo;
+
+SIGap2 = circshift(SIGap,-k3);
+SIGam2 = circshift(SIGam,k3);
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot([0:w-1]-wo,log(abs(SIGap2(wo+pp,:))),'k--','LineWidth',3,'MarkerSize',6)
+plot([0:w-1]-wo,log(abs(fDIp(wo+pp,:))),'mo-','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/SIMexpt/WoFilterCenter.m b/SIMexpt/WoFilterCenter.m
new file mode 100644
index 0000000..8299473
--- /dev/null
+++ b/SIMexpt/WoFilterCenter.m
@@ -0,0 +1,57 @@
+function [FiSMaof,NoisePower] = WoFilterCenter(FiSMao,OTFo,co,OBJpara,SFo)
+% Aim: Wiener Filtering the central frequency component
+% INPUT VARIABLES
+% FiSMao: noisy central frequency component
+% OTFo: system OTF
+% co: Wiener filter constant [=1, for minimum RMS estimate]
+% SFo: scaling factor (not significant here, so set to 1)
+% OUTPUT VARIABLES
+% FiSMaof: Wiener Filtered estimate of FiSMao
+% NoisePower: avg. noise power in FiSMao
+
+w = size(FiSMao,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Ro = sqrt( (X-wo).^2 + (Y-wo).^2 );
+
+OTFpower = OTFo.*conj(OTFo);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = FiSMao.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Object Power determination
+Aobj = OBJpara(1);
+Bobj = OBJpara(2);
+Ro(wo+1,wo+1) = 1;
+OBJpower = Aobj*(Ro.^Bobj);
+OBJpower = OBJpower.^2;
+
+%% Wiener Filtering
+FiSMaof = FiSMao.*(SFo.*conj(OTFo)./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+
+%% for cross-checking filtered estimate visually
+%{
+WFilter = (SFo.^2.*OTFpower./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+FiSMao1 = FiSMaof.*OTFo;
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(FiSMaof(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao(wo+pp,:))),'o--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao1(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,10*WFilter(wo+pp,:),'k*--','LineWidth',2,'MarkerSize',4)
+title('WoFilter')
+legend('FiSMaof','FiSMao','FiSMao1')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/SIMexpt/WoFilterSideLobe.m b/SIMexpt/WoFilterSideLobe.m
new file mode 100644
index 0000000..fb717c5
--- /dev/null
+++ b/SIMexpt/WoFilterSideLobe.m
@@ -0,0 +1,54 @@
+function [FiSMaof,NoisePower] = WoFilterSideLobe(FiSMao,OTFo,co,OBJsideP,SFo)
+% Aim: Wiener Filtering the off-center frequency component
+% INPUT VARIABLES
+% FiSMao: noisy off-center frequency component
+% OTFo: system OTF
+% co: Wiener filter constant [=1, for minimum RMS estimate]
+% OBJsideP: avg. object spectrum
+% SFo: modulation factor
+% OUTPUT VARIABLES
+% FiSMaof: Wiener Filtered estimate of FiSMao
+% NoisePower: avg. noise power in FiSMao
+
+w = size(FiSMao,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Ro = sqrt( (X-wo).^2 + (Y-wo).^2 );
+
+OTFpower = OTFo.*conj(OTFo);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = FiSMao.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Object Power determination
+OBJpower = OBJsideP.^2;
+
+%% Wiener Filtering
+FiSMaof = FiSMao.*(SFo.*conj(OTFo)./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+
+%% for cross-checking filtered estimate visually
+%{
+WFilter = (SFo.^2.*OTFpower./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+FiSMao1 = FiSMaof.*OTFo;
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(FiSMaof(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao(wo+pp,:))),'o--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao1(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,10*WFilter(wo+pp,:),'k*--','LineWidth',2,'MarkerSize',4)
+title('WoFilter')
+legend('FiSMaof','FiSMao','FiSMao1')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/SIMexpt/sim01z4.tif b/SIMexpt/sim01z4.tif
new file mode 100644
index 0000000..e8fa42d
Binary files /dev/null and b/SIMexpt/sim01z4.tif differ
diff --git a/TIRFbasic/ApodizationFunction.m b/TIRFbasic/ApodizationFunction.m
new file mode 100644
index 0000000..172f322
--- /dev/null
+++ b/TIRFbasic/ApodizationFunction.m
@@ -0,0 +1,30 @@
+function [FsumA] = ApodizationFunction(Fsum,k2fa,k2fb,k2fc,Kotf,Index)
+
+t = size(Fsum,1);
+
+OTFo_mask = OTFmaskShifted(0.*k2fa,Kotf,t);
+OTFap_mask = OTFmaskShifted(k2fa,Kotf,t);
+OTFam_mask = OTFmaskShifted(-k2fa,Kotf,t);
+OTFbp_mask = OTFmaskShifted(k2fb,Kotf,t);
+OTFbm_mask = OTFmaskShifted(-k2fb,Kotf,t);
+OTFcp_mask = OTFmaskShifted(k2fc,Kotf,t);
+OTFcm_mask = OTFmaskShifted(-k2fc,Kotf,t);
+ApoMask = OTFo_mask.*OTFap_mask.*OTFam_mask.*OTFbp_mask.*OTFbm_mask.*OTFcp_mask.*OTFcm_mask;
+clear OTFoz OTFo_mask OTFap_mask OTFam_mask OTFbp_mask OTFbm_mask OTFcp_mask OTFcm_mask
+%{
+figure;
+imshow(ApoMask,[ ])
+%}
+
+DistApoMask = bwdist(ApoMask);
+maxApoMask = max(max(DistApoMask));
+ApoFunc = double(DistApoMask./maxApoMask).^Index;
+%{
+figure;
+imshow(DistApoMask,[ ])
+figure;
+mesh(double(ApoFunc))
+%}
+
+FsumA = Fsum.*ApoFunc;
+clear ApoMask DistApoMask ApoFunc
\ No newline at end of file
diff --git a/TIRFbasic/Ifreq2opt.m b/TIRFbasic/Ifreq2opt.m
new file mode 100644
index 0000000..1b1a9ea
--- /dev/null
+++ b/TIRFbasic/Ifreq2opt.m
@@ -0,0 +1,49 @@
+function [CCop] = Ifreq2opt(k2fa,fAo,fAp,OTFo)
+% Aim: compute cross-correlation between central and off-central frequency
+% components
+% INPUT VARIABLES:
+% k2fa: nearest pixel approximation to illumination frequency vector
+% fAo: central frequency component
+% fAp: off-center frequency component
+% OTFo: system OTF
+
+w = size(OTFo,1);
+wo = w/2;
+
+% suppressing residual low frequency peaks of off-center components using
+fAo = fAo.*(1-1*OTFo.^10);
+fAp = fAp.*(1-1*OTFo.^10);
+
+fAoT = fAo.*conj(OTFo);
+fApT = fAp.*conj(OTFo);
+
+Kotf = OTFedgeF(OTFo);
+DoubleMatSize = 0;
+if ( 2*Kotf > wo )
+ DoubleMatSize = 1; % 1 for doubling fourier domain size, 0 for keeping it unchanged
+end
+
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ fAoT_temp = zeros(t,t);
+ fApT_temp = zeros(t,t);
+ fAoT_temp(wo+1:w+wo,wo+1:w+wo) = fAoT;
+ fApT_temp(wo+1:w+wo,wo+1:w+wo) = fApT;
+ clear fAoT fApT
+ fAoT = fAoT_temp;
+ fApT = fApT_temp;
+ clear fAoT_temp fApT_temp
+else
+ t = w;
+end
+to = t/2;
+u = linspace(0,t-1,t);
+v = linspace(0,t-1,t);
+[U,V] = meshgrid(u,v);
+ApT = exp( +1i.*2*pi*( k2fa(2)/t.*(U-to)+k2fa(1)/t.*(V-to) ) ).*ifft2(fApT);
+fApT0 = fft2( ApT );
+
+mA = sum(sum( fAoT.*conj(fApT0) ));
+mA = mA./sum(sum( fApT0.*conj(fApT0) ));
+
+CCop = -abs(mA);
\ No newline at end of file
diff --git a/TIRFbasic/IlluminationFreqTIRF.m b/TIRFbasic/IlluminationFreqTIRF.m
new file mode 100644
index 0000000..63ac30c
--- /dev/null
+++ b/TIRFbasic/IlluminationFreqTIRF.m
@@ -0,0 +1,71 @@
+function k2fa1 = IlluminationFreqTIRF(fDo,fDp,OTFo,k2o,thetaA)
+
+% AIM: To estimate illumination frequency vector
+% INPUT VARIABLES:
+% fDo: noisy central frequency component
+% fDp: noisy off-center frequency component
+% OTFo: system OTF
+% k2o: tentative magnitude of illumination frequency vector
+% thetaA: angular orientation of structured illumination
+
+%% for noise suppression
+fAp0 = fDp.*OTFo;
+fAo0 = fDo.*OTFo;
+
+% OTF cut-off freq
+Kotf = OTFedgeF(OTFo);
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+Zo = double( Ro<Kotf );
+
+Ka = k2o.*[cos(thetaA) sin(thetaA)];
+Ck = Ka(1) + 1i*Ka(2);
+rp = abs(Cv - Ck);
+rm = abs(Cv + Ck);
+r1 = 20; % local search radius
+Zmask = (rp>r1).*(rm>r1);
+Zmask = 1 - Zmask;
+
+% computing cross-correlation between central and off-central frequency
+% components to find nearest pixel approximation to illumination frequency
+% vector
+Kmax = k2o + r1;
+RHO = zeros(w,w);
+for p = wo-Kmax:wo+Kmax
+ for q = wo-Kmax:wo+Kmax
+ if ( Zmask(p,q) > 0 )
+ Cp = (p-wo) + 1i*(q-wo);
+ Rp = abs(Cv-Cp);
+ Zp = double( Rp<Kotf );
+ Z1 = Zo.*Zp;
+ fAp1 = circshift(fAp0,[(p-wo) (q-wo)]);
+ CC = sum(sum( fAo0.*conj(fAp1) ));
+ RHO(p,q) = CC./sum(sum(Z1));
+ end
+ end
+end
+% figure, mesh(abs(RHO))
+
+[Mx Idx] = max(abs(RHO),[],1);
+[Mx0 Idx0] = max(Mx);
+[My Idy] = max(abs(RHO),[],2);
+[My0 Idy0] = max(My);
+px0 = Idx0 - wo;
+py0 = Idy0 - wo;
+k2fa = [py0 px0]; % nearest pixel approximation
+
+%% subpixel approximation
+Ifreq2opt0 = @(k2fa0)Ifreq2opt(k2fa0,fDo,fDp,OTFo);
+% options = optimset('LargeScale','off','Algorithm',...
+% 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxIter',200,'Display','iter');
+k2fa0 = k2fa;
+[k2fa1,fval] = fminsearch(Ifreq2opt0,k2fa0,options);
+k2a = sqrt(k2fa1*k2fa1')
diff --git a/TIRFbasic/Kai2Opt.m b/TIRFbasic/Kai2Opt.m
new file mode 100644
index 0000000..319d2ae
--- /dev/null
+++ b/TIRFbasic/Kai2Opt.m
@@ -0,0 +1,38 @@
+function CCo = Kai2Opt(phaseShift,FcS1aT,FcS2aT,OTFo)
+% Aim: Computes normalized cross-power spectrum between separated frequency
+% components
+% INPUT VARIABLES:
+% phaseShift: relative phases [phaseShift2,phaseShift3]
+% FcS1aT: FT of difference between raw images 1 and 2
+% FcS2aT: FT of difference between raw images 1 and 3
+% OTFo: system OTF
+
+phaseShift2 = phaseShift(1);
+phaseShift3 = phaseShift(2);
+
+phaseShift1 = 0;
+phase2 = exp(-1i*phaseShift2) - exp(-1i*phaseShift1);
+phase3 = exp(-1i*phaseShift3) - exp(-1i*phaseShift1);
+%% Transformation Matrix
+M = [phase2 conj(phase2);
+ phase3 conj(phase3)];
+
+%% Separting the components
+%===========================================================
+Minv = inv(M);
+
+FiSMap = Minv(1,1)*FcS1aT + Minv(1,2)*FcS2aT;
+FiSMam = Minv(2,1)*FcS1aT + Minv(2,2)*FcS2aT;
+
+%{
+figure;
+mesh(log(abs(FiSMap)))
+figure;
+mesh(log(abs(FiSMam)))
+kk
+%}
+
+FiSMap0 = FiSMap.*conj(OTFo);
+FiSMam0 = FiSMam.*conj(OTFo);
+
+CCo = abs( sum(sum( FiSMap0.*conj(FiSMam0) )) );
diff --git a/TIRFbasic/MergingHeptaletsF.m b/TIRFbasic/MergingHeptaletsF.m
new file mode 100644
index 0000000..0810d21
--- /dev/null
+++ b/TIRFbasic/MergingHeptaletsF.m
@@ -0,0 +1,67 @@
+function [Fsum,Fperi,Fcent] = MergingHeptaletsF(fDIo,fDIp,fDIm,fDBo,fDBp,fDBm,...
+ fDCo,fDCp,fDCm,Ma,Mb,Mc,npDIo,npDIp,npDIm,...
+ npDBo,npDBp,npDBm,npDCo,npDCp,npDCm,k2fa,k2fb,k2fc,OBJpara,OTFo)
+
+% AIM: To merge all 9 frequency components into one using generalized
+% Wiener Filter
+% INPUT VARIABLES
+% [fDIo fDIp fDIm
+% fDBo fDBp fDBm
+% fDCo fDCp fDCm]: nine frequency components
+% Ma,Mb,Mc: modulation factors for the three illumination orientations
+% [npDIo npDIp npDIm
+% npDBo npDBp npDBm
+% npDCo,npDCp,npDCm]: noise powers corresponding to nine frequency
+% components
+% k2fa,k2fb,k2fc: illumination frequency vectors for the three
+% illumination orientations
+% OBJpara: Object spectrum parameters
+% OTFo: system OTF
+% OUTPUT VARIABLES
+% Fsum: all nine frequency components merged into one using
+% generalised Wiener Filter
+% Fperi: six off-center frequency components merged into one using
+% generalised Wiener Filter
+% Fcent: averaged of the three central frequency components
+
+% obtain signal spectrums corresponding to central and off-center
+% frequency components
+[SIGao,SIGam2,SIGap2] = TripletSNR0(OBJpara,k2fa,OTFo,fDIm);
+[SIGbo,SIGbm2,SIGbp2] = TripletSNR0(OBJpara,k2fb,OTFo,fDBm);
+[SIGco,SIGcm2,SIGcp2] = TripletSNR0(OBJpara,k2fc,OTFo,fDCm);
+
+SIGap2 = Ma*SIGap2;
+SIGam2 = Ma*SIGam2;
+SIGbp2 = Mb*SIGbp2;
+SIGbm2 = Mb*SIGbm2;
+SIGcp2 = Mc*SIGcp2;
+SIGcm2 = Mc*SIGcm2;
+
+%% Generalized Wiener-Filter computation
+SNRao = SIGao.*conj(SIGao)./npDIo;
+SNRap = SIGap2.*conj(SIGap2)./npDIp;
+SNRam = SIGam2.*conj(SIGam2)./npDIm;
+
+SNRbo = SIGbo.*conj(SIGbo)./npDBo;
+SNRbp = SIGbp2.*conj(SIGbp2)./npDBp;
+SNRbm = SIGbm2.*conj(SIGbm2)./npDBm;
+
+SNRco = SIGco.*conj(SIGco)./npDCo;
+SNRcp = SIGcp2.*conj(SIGcp2)./npDCp;
+SNRcm = SIGcm2.*conj(SIGcm2)./npDCm;
+
+ComDeno = 0.01 + ( SNRao + SNRap + SNRam + SNRao + SNRbp + SNRbm + SNRco + SNRcp + SNRcm );
+Fsum = fDIo.*SNRao + fDIp.*SNRap + fDIm.*SNRam...
+ + fDBo.*SNRbo + fDBp.*SNRbp + fDBm.*SNRbm...
+ + fDCo.*SNRco + fDCp.*SNRcp + fDCm.*SNRcm;
+Fsum = Fsum./ComDeno;
+
+ComPeri = 0.01 + ( SNRap + SNRam + SNRbp + SNRbm + SNRcp + SNRcm );
+Fperi = fDIp.*SNRap + fDIm.*SNRam + fDBp.*SNRbp + fDBm.*SNRbm...
+ + fDCp.*SNRcp + fDCm.*SNRcm;
+Fperi = Fperi./ComPeri;
+
+% averaged central frequency component
+Fcent = (fDIo + fDBo + fDCo)/3;
+
+
diff --git a/TIRFbasic/ModulationFactor.m b/TIRFbasic/ModulationFactor.m
new file mode 100644
index 0000000..2ac2dda
--- /dev/null
+++ b/TIRFbasic/ModulationFactor.m
@@ -0,0 +1,86 @@
+function Mm = ModulationFactor(fDp,k1a,OBJparaA,OTFo)
+% AIM: Determination of modulation factor
+% INPUT VARIABLES
+% fDp: off-center frequency component
+% k1a: illumination frequency vector
+% OBJparaA: Object power parameters
+% OTFo: system OTF
+% OUTPUT VARIABLE
+% Mm: modulation factor
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% magnitude of illumination vector
+k2 = sqrt(k1a*k1a');
+
+% vector along illumination direction
+kv = k1a(2) + 1i*k1a(1);
+Rp = abs(Cv+kv);
+
+% Object power parameters
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+
+% Object spectrum
+OBJp = Aobj*(Rp+0).^Bobj;
+
+% illumination vector rounded to nearest pixel
+k3 = -round(k1a);
+
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+
+% signal spectrum
+SIGap = OBJp.*OTFo;
+
+% figure;
+% mesh(log(abs(OBJp)))
+% figure;
+% mesh(log(abs(fDp)))
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = fDp.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Noise free object power computation
+Fpower = fDp.*conj(fDp) - NoisePower;
+fDp = sqrt(abs(Fpower));
+
+% frequency range over which signal power matching is done to estimate
+% modulation factor
+Zmask = (Ro > 0.2*k2).*(Ro < 0.8*k2).*(Rp > 0.2*k2);
+
+% least square approximation for modulation factor
+Mm = sum(sum(SIGap.*abs(fDp).*Zmask));
+Mm = Mm./sum(sum(SIGap.^2.*Zmask))
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(fDp(wo+pp,:))),'o-','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(SIGap(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,log(abs(Mm*SIGap(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+grid on
+box on
+kk
+%}
+
+
+
diff --git a/TIRFbasic/OBJparaOpt.m b/TIRFbasic/OBJparaOpt.m
new file mode 100644
index 0000000..8acbbbb
--- /dev/null
+++ b/TIRFbasic/OBJparaOpt.m
@@ -0,0 +1,47 @@
+function Esum = OBJparaOpt(OBJpara0,fDIoTnoisy,OTFo)
+% AIM: Determined Sum of Squared Errors (SSE) between `actual signal power' and
+% `approximated signal power'
+% INPUT VARIABLES
+% OBJpara0: [Aobj Bobj], Object power parameters
+% fDIoTnoisy: FT of central frequency component
+% OTFo: system OTF
+% OUTPUT VARIABLE
+% Esum: SSE between `actual signal spectrum' and `approximated signal spectrum'
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+Ro(wo+1,wo+1) = 1; % to avoid nan
+
+% approximated signal power calculation
+Aobj = OBJpara0(1);
+Bobj = OBJpara0(2);
+OBJpower = Aobj*(Ro.^Bobj);
+SIGpower = OBJpower.*OTFo;
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% range of frequency over which SSE is computed
+Zloop = (Ro<0.75*Kotf).*(Ro>0.25*Kotf);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = fDIoTnoisy.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Noise free object power computation
+Fpower = fDIoTnoisy.*conj(fDIoTnoisy) - NoisePower;
+fDIoTnoisy = sqrt(abs(Fpower));
+
+% SSE computation
+Error = fDIoTnoisy - SIGpower;
+Esum = sum(sum((Error.^2./Ro).*Zloop));
+
diff --git a/TIRFbasic/OBJpowerPara.m b/TIRFbasic/OBJpowerPara.m
new file mode 100644
index 0000000..9ef86ce
--- /dev/null
+++ b/TIRFbasic/OBJpowerPara.m
@@ -0,0 +1,48 @@
+function OBJparaA = OBJpowerPara(fDIoTnoisy,OTFo)
+% AIM: determination of object power parameters Aobj and Bobj
+% INPUT VARIABLES
+% fDIoTnoisy: FT of central frequency component
+% OTFo: system OTFo
+% OUTPUT VARIABLES
+% OBJparaA: [Aobj Bobj]
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+%% object power parameters through optimization
+OBJparaOpt0 = @(OBJpara0)OBJparaOpt(OBJpara0,fDIoTnoisy,OTFo);
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+
+% obtaining crude initial guesses for Aobj and Bobj
+Zm = (Ro>0.3*Kotf).*(Ro<0.4*Kotf);
+Aobj = sum(sum(abs(fDIoTnoisy.*Zm)))./sum(sum(Zm));
+Bobj = -0.5;
+OBJpara0 = [Aobj Bobj];
+
+% optimization step
+[OBJparaA,fval] = fminsearch(OBJparaOpt0,OBJpara0,options);
+
+%% plot for cross-checking the result
+%{
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+OBJo = Aobj*(Ro.^Bobj);
+SIGp = OBJo.*OTFo;
+pp = 3;
+figure;
+hold on
+plot([0:w-1]-wo,log(abs(fDIoTnoisy(wo+pp,:))),'k--','LineWidth',3,'MarkerSize',6)
+plot([0:w-1]-wo,log(abs(SIGp(wo+pp,:))),'mo-','LineWidth',2,'MarkerSize',6)
+legend('acutal signal power','avg. signal power')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/TIRFbasic/OTFdoubling.m b/TIRFbasic/OTFdoubling.m
new file mode 100644
index 0000000..58aaa6e
--- /dev/null
+++ b/TIRFbasic/OTFdoubling.m
@@ -0,0 +1,17 @@
+function OTF2 = OTFdoubling(OTFo,DoubleMatSize)
+
+%% embeds OTFo in doubled range of frequencies
+
+w = size(OTFo,1);
+wo = w/2;
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ OTF2 = zeros(2*w,2*w);
+ OTF2(wo+1:w+wo,wo+1:w+wo) = OTFo;
+ clear OTFo
+else
+ OTF2 = OTFo;
+ clear OTFo
+end
\ No newline at end of file
diff --git a/TIRFbasic/OTFedgeF.m b/TIRFbasic/OTFedgeF.m
new file mode 100644
index 0000000..e4d6a7f
--- /dev/null
+++ b/TIRFbasic/OTFedgeF.m
@@ -0,0 +1,13 @@
+function Kotf = OTFedgeF(OTFo)
+
+w = size(OTFo,1);
+wo = w/2;
+
+OTF1 = OTFo(wo+1,:);
+OTFmax = max(max(abs(OTFo)));
+OTFtruncate = 0.01;
+i = 1;
+while ( abs(OTF1(1,i))<OTFtruncate*OTFmax )
+ Kotf = wo+1-i;
+ i = i + 1;
+end
\ No newline at end of file
diff --git a/TIRFbasic/OTFmaskShifted.m b/TIRFbasic/OTFmaskShifted.m
new file mode 100644
index 0000000..a60ac98
--- /dev/null
+++ b/TIRFbasic/OTFmaskShifted.m
@@ -0,0 +1,13 @@
+function [OTFap_mask] = OTFmaskShifted(k2fa,Kotf,w)
+% OTFap_mask = complement of shifted OTF
+
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+Iy = k2fa(1,1) + (wo+1);
+Ix = k2fa(1,2) + (wo+1);
+
+Ro = sqrt( (X+1-Ix).^2 + (Y+1-Iy).^2 );
+OTFap_mask = Ro > Kotf;
\ No newline at end of file
diff --git a/TIRFbasic/PCMfilteringF.m b/TIRFbasic/PCMfilteringF.m
new file mode 100644
index 0000000..0c68963
--- /dev/null
+++ b/TIRFbasic/PCMfilteringF.m
@@ -0,0 +1,146 @@
+function [fDof,fDp2,fDm2,npDo,npDp,npDm,Mm,DoubleMatSize]...
+ = PCMfilteringF(fDo,fDp,fDm,OTFo,OBJparaA,kA)
+
+% AIM: obtaining Wiener Filtered estimates of noisy frequency components
+% INPUT VARIABLES
+% fDo,fDp,fDm: noisy estimates of separated frequency components
+% OTFo: system OTF
+% OBJparaA: object power parameters
+% kA: illumination vector
+% OUTPUT VARIABLES
+% fDof,fDp2,fDm2: Wiener Filtered estimates of fDo,fDp,fDm
+% npDo,npDp,npDm: avg. noise power in fDo,fDp,fDm
+% Mm: illumination modulation factor
+% DoubleMatSize: parameter for doubling FT size if necessary
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% object power parameters
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+
+% Wiener Filtering central frequency component
+SFo = 1;
+co = 1.0;
+[fDof,npDo] = WoFilterCenter(fDo,OTFo,co,OBJparaA,SFo);
+
+% modulation factor determination
+Mm = ModulationFactor(fDp,kA,OBJparaA,OTFo);
+
+%% Duplex power (default)
+kv = kA(2) + 1i*kA(1); % vector along illumination direction
+Rp = abs(Cv-kv);
+Rm = abs(Cv+kv);
+OBJp = Aobj*(Rp.^Bobj);
+OBJm = Aobj*(Rm.^Bobj);
+k3 = round(kA);
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+OBJm(wo+1-k3(1),wo+1-k3(2)) = 0.25*OBJm(wo+2-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+2-k3(2))...
+ + 0.25*OBJm(wo+0-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+0-k3(2));
+
+% Filtering side lobes (off-center frequency components)
+SFo = Mm;
+[fDpf,npDp] = WoFilterSideLobe(fDp,OTFo,co,OBJm,SFo);
+[fDmf,npDm] = WoFilterSideLobe(fDm,OTFo,co,OBJp,SFo);
+
+%% doubling Fourier domain size if necessary
+DoubleMatSize = 0;
+if ( 2*Kotf > wo )
+ DoubleMatSize = 1; % 1 for doubling fourier domain size, 0 for keeping it unchanged
+end
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ to = t/2;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ [U,V] = meshgrid(u,v);
+ fDoTemp = zeros(2*w,2*w);
+ fDpTemp = zeros(2*w,2*w);
+ fDmTemp = zeros(2*w,2*w);
+ OTFtemp = zeros(2*w,2*w);
+ fDoTemp(wo+1:w+wo,wo+1:w+wo) = fDof;
+ fDpTemp(wo+1:w+wo,wo+1:w+wo) = fDpf;
+ fDmTemp(wo+1:w+wo,wo+1:w+wo) = fDmf;
+ OTFtemp(wo+1:w+wo,wo+1:w+wo) = OTFo;
+ clear fDof fDpf fDmf OTFo w wo x y X Y
+ fDof = fDoTemp;
+ fDpf = fDpTemp;
+ fDmf = fDmTemp;
+ OTFo = OTFtemp;
+ clear fDoTemp fDpTemp fDmTemp OTFtemp
+else
+ t = w;
+ to = t/2;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ [U,V] = meshgrid(u,v);
+end
+
+% Shifting the off-center frequency components to their correct location
+fDp1 = fft2(ifft2(fDpf).*exp( +1i.*2*pi*(kA(2)/t.*(U-to) + kA(1)/t.*(V-to)) ));
+fDm1 = fft2(ifft2(fDmf).*exp( -1i.*2*pi*(kA(2)/t.*(U-to) + kA(1)/t.*(V-to)) ));
+%{
+figure;
+mesh(log(abs(fDpf)))
+figure;
+mesh(log(abs(fDp1)))
+kk
+%}
+
+%% Phase matching
+Cv = (U-to) + 1i*(V-to);
+Ro = abs(Cv);
+Rp = abs(Cv-kv);
+k2 = sqrt(kA*kA');
+
+% frequency range over which corrective phase is determined
+Zmask = (Ro < 0.8*k2).*(Rp < 0.8*k2);
+
+% corrective phase
+Angle0 = angle( sum(sum( fDof.*conj(fDp1).*Zmask )) );
+
+% phase correction
+fDp2 = exp(+1i*Angle0).*fDp1;
+fDm2 = exp(-1i*Angle0).*fDm1;
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot(u-to,angle(fDof(to+pp,:)).*180/pi,'o-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(fDp1(to+pp,:)).*180/pi,'ro-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(fDp2(to+pp,:)).*180/pi,'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+figure;
+hold on
+plot3(u-to,real(fDof(to+1,:)),imag(fDof(to+1,:)),'o-','LineWidth',2,'MarkerSize',8)
+plot3(u-to,real(fDp1(to+1,:)),imag(fDp1(to+1,:)),'ro-','LineWidth',2,'MarkerSize',8)
+plot3(u-to,real(fDp2(to+1,:)),imag(fDp2(to+1,:)),'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+
+ff1 = fDof.*conj(fDp1);
+ff2 = fDof.*conj(fDp2);
+figure;
+hold on
+plot(u-to,angle(ff1(to+pp,:)).*180/pi,'ro-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(ff2(to+pp,:)).*180/pi,'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+%}
+
diff --git a/TIRFbasic/PCMseparateF.m b/TIRFbasic/PCMseparateF.m
new file mode 100644
index 0000000..70e2d56
--- /dev/null
+++ b/TIRFbasic/PCMseparateF.m
@@ -0,0 +1,57 @@
+function [fDo,fDp,fDm]...
+ = PCMseparateF(S1aTnoisy,S2aTnoisy,S3aTnoisy,OTFo)
+% AIM: obtaining the noisy estimates of three frequency components
+% INPUT VARIABLES
+% S1aTnoisy,S2aTnoisy,S3aTnoisy: 3 raw SIM images with identical
+% illumination pattern orientation
+% but different phase shifts
+% OTFo: system OTF
+% OUTPUT VARIABLES
+% fDo,fDp,fDm: noisy estimates of separated frequency components
+
+
+w = size(S1aTnoisy,1);
+wo = w/2;
+
+% computing PSFe for edge tapering SIM images
+PSFd = real(fftshift( ifft2(fftshift(OTFo.^3)) ));
+PSFd = PSFd/max(max(PSFd));
+PSFd = PSFd/sum(sum(PSFd));
+h = 30;
+PSFe = PSFd(wo-h+1:wo+h,wo-h+1:wo+h);
+
+% edge tapering raw SIM images
+S1aTnoisy = edgetaper(S1aTnoisy,PSFe);
+S2aTnoisy = edgetaper(S2aTnoisy,PSFe);
+S3aTnoisy = edgetaper(S3aTnoisy,PSFe);
+
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+fS2aTnoisy = fftshift(fft2(S2aTnoisy));
+fS3aTnoisy = fftshift(fft2(S3aTnoisy));
+
+fDuplex2 = fS2aTnoisy - fS1aTnoisy;
+fDuplex3 = fS3aTnoisy - fS1aTnoisy;
+
+%% Optimizing the relative phases
+Kai2Opt0 = @(phase0)Kai2Opt(phase0,fDuplex2,fDuplex3,OTFo);
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+phase0 = [pi/2 -pi/2]; % initial guess
+[phaseShift,fval] = fminsearch(Kai2Opt0,phase0,options);
+phaseShift*180/pi
+
+%% Separating the three frequency components
+phaseShift0 = 0;
+[fDo,fDp,fDm] = SeparatedComponents2D(...
+ phaseShift,phaseShift0,fS1aTnoisy,fS2aTnoisy,fS3aTnoisy);
+
+%% for visual verification
+%{
+figure;
+mesh(log(abs(fDo)))
+figure;
+mesh(log(abs(fDp)))
+figure;
+mesh(log(abs(fDm)))
+% kk
+%}
\ No newline at end of file
diff --git a/TIRFbasic/PsfOtf.m b/TIRFbasic/PsfOtf.m
new file mode 100644
index 0000000..0120d1c
--- /dev/null
+++ b/TIRFbasic/PsfOtf.m
@@ -0,0 +1,34 @@
+function [yy0,OTF2dc] = PsfOtf(w,scale)
+
+% AIM: To generate PSF and OTF using Bessel function
+% INPUT VARIABLES
+% w: image size
+% scale: a parameter used to adjust PSF/OTF width
+% OUTPUT VRAIBLES
+% yyo: system PSF
+% OTF2dc: system OTF
+
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+%% Generation of the PSF with Besselj.
+% scale = 0.3;
+R=sqrt(min(X,abs(X-w)).^2+min(Y,abs(Y-w)).^2);
+yy=abs(2*besselj(1,scale*R+eps,1)./(scale*R+eps)).^2; %0.5 is introduced to make PSF wider
+yy0 = fftshift(yy);
+% figure;
+% mesh(X,Y,yy0)
+% figure;
+% imshow(yy0,[])
+% axis square
+
+%Generate 2D OTF.
+OTF2d=fft2(yy);
+OTF2dmax = max(max(abs(OTF2d)));
+OTF2d = OTF2d./OTF2dmax;
+OTF2dc = abs(fftshift(OTF2d));
+% figure;
+% mesh(abs(OTF2dc))
+% title('OTF of the system')
+% axis square
\ No newline at end of file
diff --git a/TIRFbasic/SIMimagesF.m b/TIRFbasic/SIMimagesF.m
new file mode 100644
index 0000000..dc977bf
--- /dev/null
+++ b/TIRFbasic/SIMimagesF.m
@@ -0,0 +1,281 @@
+function [S1aTnoisy S2aTnoisy S3aTnoisy ...
+ S1bTnoisy S2bTnoisy S3bTnoisy ...
+ S1cTnoisy S2cTnoisy S3cTnoisy ...
+ DIoTnoisy DIoT] = SIMimagesF(k2,...
+ DIo,PSFo,OTFo,ModFac,NoiseLevel,UsePSF)
+
+% AIM: to generate raw sim images
+% INPUT VARIABLES
+% k2: illumination frequency
+% DIo: specimen image
+% PSFo: system PSF
+% OTFo: system OTF
+% UsePSF: 1 (to blur SIM images by convloving with PSF)
+% 0 (to blur SIM images by truncating its fourier content beyond OTF)
+% NoiseLevel: percentage noise level for generating gaussian noise
+% OUTPUT VARIABLES
+% [S1aTnoisy S2aTnoisy S3aTnoisy
+% S1bTnoisy S2bTnoisy S3bTnoisy
+% S1cTnoisy S2cTnoisy S3cTnoisy]: nine raw sim images
+% DIoTnoisy: noisy wide field image
+% DIoT: noise-free wide field image
+
+w = size(DIo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+%% illunination phase shifts along the three directions
+p0Ao = 0*pi/3;
+p0Ap = 2*pi/3;
+p0Am = 4*pi/3;
+p0Bo = 0*pi/3;
+p0Bp = 2*pi/3;
+p0Bm = 4*pi/3;
+p0Co = 0*pi/3;
+p0Cp = 2*pi/3;
+p0Cm = 4*pi/3;
+
+%% Illuminating patterns
+alpha = 0*pi/6;
+% orientation direction of illumination patterns
+thetaA = 0*pi/3 + alpha;
+thetaB = 1*pi/3 + alpha;
+thetaC = 2*pi/3 + alpha;
+% illumination frequency vectors
+k2a = (k2/w).*[cos(thetaA) sin(thetaA)];
+k2b = (k2/w).*[cos(thetaB) sin(thetaB)];
+k2c = (k2/w).*[cos(thetaC) sin(thetaC)];
+% -------------------------------------------------------
+% mean illumination intensity
+mA = 0.5;
+mB = 0.5;
+mC = 0.5;
+% amplitude of illumination intensity above mean
+aA = 0.5*ModFac;
+aB = 0.5*ModFac;
+aC = 0.5*ModFac;
+
+% random phase shift errors
+NN = 2*(0.5-rand(9,1))*pi/18;
+
+% illunination phase shifts with random errors
+psAo = p0Ao + NN(1,1);
+psAp = p0Ap + NN(2,1);
+psAm = p0Am + NN(3,1);
+psBo = p0Bo + NN(4,1);
+psBp = p0Bp + NN(5,1);
+psBm = p0Bm + NN(6,1);
+psCo = p0Co + NN(7,1);
+psCp = p0Cp + NN(8,1);
+psCm = p0Cm + NN(9,1);
+[psAo psAp psAm;
+ psBo psBp psBm;
+ psCo psCp psCm]*180/pi
+
+% illunination patterns
+sAo = mA + aA*cos(2*pi*(k2a(1,1).*(X-wo)+k2a(1,2).*(Y-wo))+psAo); % illuminated signal (0 phase)
+sAp = mA + aA*cos(2*pi*(k2a(1,1).*(X-wo)+k2a(1,2).*(Y-wo))+psAp); % illuminated signal (+ phase)
+sAm = mA + aA*cos(2*pi*(k2a(1,1).*(X-wo)+k2a(1,2).*(Y-wo))+psAm); % illuminated signal (- phase)
+sBo = mB + aB*cos(2*pi*(k2b(1,1).*(X-wo)+k2b(1,2).*(Y-wo))+psBo); % illuminated signal (0 phase)
+sBp = mB + aB*cos(2*pi*(k2b(1,1).*(X-wo)+k2b(1,2).*(Y-wo))+psBp); % illuminated signal (+ phase)
+sBm = mB + aB*cos(2*pi*(k2b(1,1).*(X-wo)+k2b(1,2).*(Y-wo))+psBm); % illuminated signal (- phase)
+sCo = mC + aC*cos(2*pi*(k2c(1,1).*(X-wo)+k2c(1,2).*(Y-wo))+psCo); % illuminated signal (0 phase)
+sCp = mC + aC*cos(2*pi*(k2c(1,1).*(X-wo)+k2c(1,2).*(Y-wo))+psCp); % illuminated signal (+ phase)
+sCm = mC + aC*cos(2*pi*(k2c(1,1).*(X-wo)+k2c(1,2).*(Y-wo))+psCm); % illuminated signal (- phase)
+
+%{
+fftemp = fftshift(fft2(sAo));
+figure;
+hold on
+plot(x-wo,abs(OTFo(wo+1,:)),'--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,0.5*abs(fftemp(wo+1,:))./max(max(abs(fftemp))),'r--','LineWidth',2,'MarkerSize',6)
+legend('OTFo','illumination spectrum')
+grid on
+box on
+%}
+
+%% superposed Objects
+s1a = DIo.*sAo; % superposed signal (0 phase)
+s2a = DIo.*sAp; % superposed signal (+ phase)
+s3a = DIo.*sAm; % superposed signal (- phase)
+s1b = DIo.*sBo;
+s2b = DIo.*sBp;
+s3b = DIo.*sBm;
+s1c = DIo.*sCo;
+s2c = DIo.*sCp;
+s3c = DIo.*sCm;
+%{
+figure;
+subplot(3,3,1)
+imshow(s1a,[ ])
+title('s1a')
+subplot(3,3,4)
+imshow(s2a,[ ])
+title('s2a')
+subplot(3,3,7)
+imshow(s3a,[ ])
+title('s3a')
+subplot(3,3,2)
+imshow(s1b,[ ])
+title('s1b')
+subplot(3,3,5)
+imshow(s2b,[ ])
+title('s2b')
+subplot(3,3,8)
+imshow(s3b,[ ])
+title('s3b')
+subplot(3,3,3)
+imshow(s1c,[ ])
+title('s1c')
+subplot(3,3,6)
+imshow(s2c,[ ])
+title('s2c')
+subplot(3,3,9)
+imshow(s3c,[ ])
+title('s3c')
+%}
+
+
+%% superposed (noise-free) Images
+PSFsum = sum(sum(PSFo));
+if ( UsePSF == 1 )
+ DIoT = conv2(DIo,PSFo,'same')./PSFsum;
+ S1aT = conv2(s1a,PSFo,'same')./PSFsum;
+ S2aT = conv2(s2a,PSFo,'same')./PSFsum;
+ S3aT = conv2(s3a,PSFo,'same')./PSFsum;
+ S1bT = conv2(s1b,PSFo,'same')./PSFsum;
+ S2bT = conv2(s2b,PSFo,'same')./PSFsum;
+ S3bT = conv2(s3b,PSFo,'same')./PSFsum;
+ S1cT = conv2(s1c,PSFo,'same')./PSFsum;
+ S2cT = conv2(s2c,PSFo,'same')./PSFsum;
+ S3cT = conv2(s3c,PSFo,'same')./PSFsum;
+else
+ DIoT = ifft2( fft2(DIo).*fftshift(OTFo) );
+ S1aT = ifft2( fft2(s1a).*fftshift(OTFo) );
+ S2aT = ifft2( fft2(s2a).*fftshift(OTFo) );
+ S3aT = ifft2( fft2(s3a).*fftshift(OTFo) );
+ S1bT = ifft2( fft2(s1b).*fftshift(OTFo) );
+ S2bT = ifft2( fft2(s2b).*fftshift(OTFo) );
+ S3bT = ifft2( fft2(s3b).*fftshift(OTFo) );
+ S1cT = ifft2( fft2(s1c).*fftshift(OTFo) );
+ S2cT = ifft2( fft2(s2c).*fftshift(OTFo) );
+ S3cT = ifft2( fft2(s3c).*fftshift(OTFo) );
+
+ DIoT = real(DIoT);
+ S1aT = real(S1aT);
+ S2aT = real(S2aT);
+ S3aT = real(S3aT);
+ S1bT = real(S1bT);
+ S2bT = real(S2bT);
+ S3bT = real(S3bT);
+ S1cT = real(S1cT);
+ S2cT = real(S2cT);
+ S3cT = real(S3cT);
+end
+%{
+figure;
+subplot(3,3,1)
+imshow(S1aT,[ ])
+title('S1aT')
+subplot(3,3,4)
+imshow(S2aT,[ ])
+title('S2aT')
+subplot(3,3,7)
+imshow(S3aT,[ ])
+title('S3aT')
+subplot(3,3,2)
+imshow(S1bT,[ ])
+title('S1bT')
+subplot(3,3,5)
+imshow(S2bT,[ ])
+title('S2bT')
+subplot(3,3,8)
+imshow(S3bT,[ ])
+title('S3bT')
+subplot(3,3,3)
+imshow(S1cT,[ ])
+title('S1cT')
+subplot(3,3,6)
+imshow(S2cT,[ ])
+title('S2cT')
+subplot(3,3,9)
+imshow(S3cT,[ ])
+title('S3cT')
+%}
+
+
+%% Gaussian noise generation
+aNoise = NoiseLevel/100; % corresponds to 10% noise
+% SNR = 1/aNoise
+% SNRdb = 20*log10(1/aNoise)
+nDIoT = random('norm', 0, aNoise*std2(DIoT), w , w);
+nS1aT = random('norm', 0, aNoise*std2(S1aT), w , w);
+nS2aT = random('norm', 0, aNoise*std2(S2aT), w , w);
+nS3aT = random('norm', 0, aNoise*std2(S3aT), w , w);
+nS1bT = random('norm', 0, aNoise*std2(S2bT), w , w);
+nS2bT = random('norm', 0, aNoise*std2(S2bT), w , w);
+nS3bT = random('norm', 0, aNoise*std2(S2bT), w , w);
+nS1cT = random('norm', 0, aNoise*std2(S3cT), w , w);
+nS2cT = random('norm', 0, aNoise*std2(S3cT), w , w);
+nS3cT = random('norm', 0, aNoise*std2(S3cT), w , w);
+
+%% noise added raw SIM images
+NoiseFrac = 1; %may be set to 0 to avoid noise addition
+DIoTnoisy = DIoT + NoiseFrac*nDIoT;
+S1aTnoisy = S1aT + NoiseFrac*nS1aT;
+S2aTnoisy = S2aT + NoiseFrac*nS2aT;
+S3aTnoisy = S3aT + NoiseFrac*nS3aT;
+S1bTnoisy = S1bT + NoiseFrac*nS1bT;
+S2bTnoisy = S2bT + NoiseFrac*nS2bT;
+S3bTnoisy = S3bT + NoiseFrac*nS3bT;
+S1cTnoisy = S1cT + NoiseFrac*nS1cT;
+S2cTnoisy = S2cT + NoiseFrac*nS2cT;
+S3cTnoisy = S3cT + NoiseFrac*nS3cT;
+%{
+figure;
+subplot(3,3,1)
+imshow(S1aTnoisy,[ ])
+title('S1aTnoisy')
+subplot(3,3,4)
+imshow(S2aTnoisy,[ ])
+title('S2aTnoisy')
+subplot(3,3,7)
+imshow(S3aTnoisy,[ ])
+title('S3aTnoisy')
+subplot(3,3,2)
+imshow(S1bTnoisy,[ ])
+title('S1bTnoisy')
+subplot(3,3,5)
+imshow(S2bTnoisy,[ ])
+title('S2bTnoisy')
+subplot(3,3,8)
+imshow(S3bTnoisy,[ ])
+title('S3bTnoisy')
+subplot(3,3,3)
+imshow(S1cTnoisy,[ ])
+title('S1cTnoisy')
+subplot(3,3,6)
+imshow(S2cTnoisy,[ ])
+title('S2cTnoisy')
+subplot(3,3,9)
+imshow(S3cTnoisy,[ ])
+title('S3cTnoisy')
+%}
+
+%{
+figure;
+imshow(S1aT,[ ])
+title('S1aT')
+figure;
+imshow(S1aTnoisy,[ ])
+title('S1aTnoisy')
+
+figure;
+hold on
+plot(x-wo,S1aT(wo+1,:),'o-','LineWidth',2,'MarkerSize',6)
+plot(x-wo,S1aTnoisy(wo+1,:),'r*--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/TIRFbasic/SIMo.m b/TIRFbasic/SIMo.m
new file mode 100644
index 0000000..fbcc4c8
--- /dev/null
+++ b/TIRFbasic/SIMo.m
@@ -0,0 +1,104 @@
+clear all
+close all
+clc
+
+w = 512;
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+%% Generation of the PSF with Besselj.
+scale = 0.63; % used to adjust PSF/OTF width
+[PSFo,OTFo] = PsfOtf(w,scale);
+%{
+figure;
+mesh(PSFo)
+title('system PSF')
+
+figure;
+mesh(OTFo)
+title('system OTF')
+%}
+
+%% Reading input file
+Io1 = imread('testpat.tiff');
+Io = Io1(257:768,257:768); % selecting central 512x512 of image
+DIo = double(Io);
+figure;
+imshow(DIo,[])
+
+
+%% Generating raw SIM Images
+k2 = 120.23; % illumination freq
+ModFac = 0.8; % modulation factor
+NoiseLevel = 10; % in percentage
+UsePSF = 0; % 1(to blur using PSF), 0(to blur using OTF)
+[S1aTnoisy S2aTnoisy S3aTnoisy ...
+ S1bTnoisy S2bTnoisy S3bTnoisy ...
+ S1cTnoisy S2cTnoisy S3cTnoisy ...
+ DIoTnoisy DIoT] = SIMimagesF(k2,DIo,PSFo,OTFo,ModFac,NoiseLevel,UsePSF);
+%{
+figure;
+imshow(S1aTnoisy,[])
+title('Raw SIM image')
+figure;
+imshow(S1bTnoisy,[])
+title('Raw SIM image')
+figure;
+imshow(S1cTnoisy,[])
+title('Raw SIM image')
+figure;
+imshow(DIoTnoisy,[])
+title('Noisy Wide-Field image')
+figure;
+imshow(DIoT,[])
+title('Noise-free Wide-Field image')
+break
+%}
+
+% tentative values to define search region for illumination frequency
+% determination
+k2o = 125; % tentative illumination frequency
+thetaA = 0*pi/3; % orientations of structured illumination
+thetaB = 1*pi/3;
+thetaC = 2*pi/3;
+
+%% obtaining the noisy estimates of three frequency components
+[fAo,fAp,fAm]...
+ = PCMseparateF(S1aTnoisy,S2aTnoisy,S3aTnoisy,OTFo);
+kA = IlluminationFreqTIRF(fAo,fAp,OTFo,k2o,thetaA);
+[fBo,fBp,fBm]...
+ = PCMseparateF(S1bTnoisy,S2bTnoisy,S3bTnoisy,OTFo);
+kB = IlluminationFreqTIRF(fBo,fBp,OTFo,k2o,thetaB);
+[fCo,fCp,fCm]...
+ = PCMseparateF(S1cTnoisy,S2cTnoisy,S3cTnoisy,OTFo);
+kC = IlluminationFreqTIRF(fCo,fCp,OTFo,k2o,thetaC);
+
+% averaging the central frequency components
+fCent = (fAo + fBo + fCo)/3;
+
+% Object power parameters determination
+OBJparaA = OBJpowerPara(fCent,OTFo);
+
+% break
+%% Wiener Filtering the noisy frequency components
+[fAof,fApf,fAmf,Nao,Nap,Nam,Ma,DoubleMatSize]...
+ = PCMfilteringF(fAo,fAp,fAm,OTFo,OBJparaA,kA);
+[fBof,fBpf,fBmf,Nbo,Nbp,Nbm,Mb,~]...
+ = PCMfilteringF(fAo,fBp,fBm,OTFo,OBJparaA,kB);
+[fCof,fCpf,fCmf,Nco,Ncp,Ncm,Mc,~]...
+ = PCMfilteringF(fAo,fCp,fCm,OTFo,OBJparaA,kC);
+
+%% doubling Fourier domain size if necessary
+OTFo = OTFdoubling(OTFo,DoubleMatSize);
+
+%% merging all 9 frequency components using generalized Wiener Filter
+[Fsum,Fperi,Fcent] = MergingHeptaletsF(fAof,fApf,fAmf,...
+ fBof,fBpf,fBmf,fCof,fCpf,fCmf,...
+ Ma,Mb,Mc,Nao,Nap,Nam,Nbo,Nbp,Nbm,...
+ Nco,Ncp,Ncm,kA,kB,kC,OBJparaA,OTFo);
+
+% Plotting SIM results
+SIMplot(Fsum,Fperi,Fcent,OTFo,kA,kB,kC,S1aTnoisy);
+
diff --git a/TIRFbasic/SIMplot.m b/TIRFbasic/SIMplot.m
new file mode 100644
index 0000000..b83fd74
--- /dev/null
+++ b/TIRFbasic/SIMplot.m
@@ -0,0 +1,84 @@
+function [] = SIMplot(Fsum,Fperi,Fcent,OTFo,kA,kB,kC,S1aTnoisy)
+
+%% recontructed SIM images
+Dcent = real( ifft2(fftshift(Fcent)) );
+Dsum = real( ifft2(fftshift(Fsum)) );
+Dperi = real( ifft2(fftshift(Fperi)) );
+
+t = size(OTFo,1);
+h = 1*30; % for removing the image edges
+figure;
+imshow(Dcent(h+1:t-h,h+1:t-h),[])
+title('Wiener-Filtered wide-field')
+figure;
+imshow(Dsum(h+1:t-h,h+1:t-h),[])
+title('SIM image')
+figure;
+imshow(Dperi(h+1:t-h,h+1:t-h),[])
+title('SIM image (using only off-center frequency components)')
+
+% appodizing the merged frequency components
+Index = 0.4;
+Kotf = OTFedgeF(OTFo);
+[FsumA] = ApodizationFunction(Fsum,kA,kB,kC,Kotf,Index);
+[FperiA] = ApodizationFunction(Fperi,kA,kB,kC,Kotf,Index);
+DsumA = real( ifft2(fftshift(FsumA)) );
+DperiA = real( ifft2(fftshift(FperiA)) );
+
+figure;
+imshow(DsumA(h+1:t-h,h+1:t-h),[])
+title('appodized SIM image')
+figure;
+imshow(DperiA(h+1:t-h,h+1:t-h),[])
+title('appodized SIM image (using only off-center frequency components)')
+
+
+%% Frequency Plots
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+w = size(OTFo,1);
+wo = w/2;
+w1 = size(fS1aTnoisy,1);
+if ( w > w1 )
+ fS1aTnoisy0 = zeros(w,w);
+ fS1aTnoisy0(wo-w1/2+1:wo+w1/2,wo-w1/2+1:wo+w1/2) = fS1aTnoisy;
+ clear fS1aTnoisy
+ fS1aTnoisy = fS1aTnoisy0;
+ clear fS1aTnoisy0;
+end
+
+
+p = 10;
+minL1 = min(min( abs(fS1aTnoisy).^(1/p) ));
+minL2 = min(min( abs(Fcent).^(1/p) ));
+minL3 = min(min( abs(Fsum).^(1/p) ));
+minL4 = min(min( abs(Fperi).^(1/p) ));
+minL5 = min(min( abs(FsumA).^(1/p) ));
+maxL1 = max(max( abs(fS1aTnoisy).^(1/p) ));
+maxL2 = max(max( abs(Fcent).^(1/p) ));
+maxL3 = max(max( abs(Fsum).^(1/p) ));
+maxL4 = max(max( abs(Fperi).^(1/p) ));
+maxL5 = max(max( abs(FsumA).^(1/p) ));
+minL = min([minL1,minL2,minL3,minL4,minL5]);
+maxL = max([maxL1,maxL2,maxL3,maxL4,maxL5]);
+
+figure;
+imshow(abs(fS1aTnoisy).^(1/p),[minL maxL])
+title('fS1aTnoisy')
+
+figure;
+imshow(abs(Fcent).^(1/p),[minL maxL])
+title('Weiner Filtered frequency')
+figure;
+imshow(abs(Fsum).^(1/p),[minL maxL])
+title('SIM frequency')
+figure;
+imshow(abs(Fperi).^(1/p),[minL maxL])
+title('SIM (off-center frequency components)')
+figure;
+imshow(abs(FsumA).^(1/p),[minL maxL])
+title('appodized SIM frequency')
+figure;
+imshow(abs(FperiA).^(1/p),[minL maxL])
+title('appodized SIM (off-center frequency components)')
+
+
diff --git a/TIRFbasic/SeparatedComponents2D.m b/TIRFbasic/SeparatedComponents2D.m
new file mode 100644
index 0000000..d124752
--- /dev/null
+++ b/TIRFbasic/SeparatedComponents2D.m
@@ -0,0 +1,22 @@
+function [FiSMao,FiSMap,FiSMam] = SeparatedComponents2D(...
+ phaseShift,phaseShift0,FcS1aT,FcS2aT,FcS3aT)
+% Aim: Unmixing the frequency components of raw SIM images
+% phaseShift,phaseShift0: illumination phase shifts
+% FcS1aT,FcS2aT,FcS3aT: FT of raw SIM images
+% FiSMao,FiSMap,FiSMam: unmixed frequency components of raw SIM images
+
+phaseShift1 = phaseShift(1);
+phaseShift2 = phaseShift(2);
+MF = 1.0;
+%% Transformation Matrix
+M = 0.5*[1 0.5*MF*exp(-1i*phaseShift0) 0.5*MF*exp(+1i*phaseShift0);
+ 1 0.5*MF*exp(-1i*phaseShift1) 0.5*MF*exp(+1i*phaseShift1);
+ 1 0.5*MF*exp(-1i*phaseShift2) 0.5*MF*exp(+1i*phaseShift2)];
+
+%% Separting the components
+%===========================================================
+Minv = inv(M);
+
+FiSMao = Minv(1,1)*FcS1aT + Minv(1,2)*FcS2aT + Minv(1,3)*FcS3aT;
+FiSMap = Minv(2,1)*FcS1aT + Minv(2,2)*FcS2aT + Minv(2,3)*FcS3aT;
+FiSMam = Minv(3,1)*FcS1aT + Minv(3,2)*FcS2aT + Minv(3,3)*FcS3aT;
diff --git a/TIRFbasic/TripletSNR0.m b/TIRFbasic/TripletSNR0.m
new file mode 100644
index 0000000..6bf5b1d
--- /dev/null
+++ b/TIRFbasic/TripletSNR0.m
@@ -0,0 +1,66 @@
+function [SIGao,SIGap2,SIGam2] = TripletSNR0(OBJpara,k2fa,OTFo,fDIp)
+
+% AIM: To obtain signal spectrums corresponding to central and off-center
+% frequency components
+% INPUT VARIABLES
+% OBJpara: object power parameters
+% k2fa: illumination frequency vector
+% OTFo: system OTF
+% fDIp: one of the off-center frequency component (utilized here only for
+% visual verification of computation)
+% OUTPUT VARIABLES
+% SIGao: signal spectrum corresponding to central frequency component
+% SIGap2,SIGam2: signal spectrums corresponding to off-center frequency components
+
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% object spectrum parameters
+Aobj = OBJpara(1);
+Bobj = OBJpara(2);
+
+% object spectrums (default)
+kv = k2fa(2) + 1i*k2fa(1); % vector along illumination direction
+Rp = abs(Cv-kv);
+Rm = abs(Cv+kv);
+OBJo = Aobj*(Ro.^Bobj);
+OBJp = Aobj*(Rp.^Bobj);
+OBJm = Aobj*(Rm.^Bobj);
+
+OBJo(wo+1,wo+1) = 0.25*OBJo(wo+2,wo+1) + 0.25*OBJo(wo+1,wo+2)...
+ + 0.25*OBJo(wo+0,wo+1) + 0.25*OBJo(wo+1,wo+0);
+
+k3 = round(k2fa);
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+OBJm(wo+1-k3(1),wo+1-k3(2)) = 0.25*OBJm(wo+2-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+2-k3(2))...
+ + 0.25*OBJm(wo+0-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+0-k3(2));
+
+% signal spectrums
+SIGao = OBJo.*OTFo;
+SIGap = OBJp.*OTFo;
+SIGam = OBJm.*OTFo;
+
+SIGap2 = circshift(SIGap,-k3);
+SIGam2 = circshift(SIGam,k3);
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot([0:w-1]-wo,log(abs(SIGap2(wo+pp,:))),'k--','LineWidth',3,'MarkerSize',6)
+plot([0:w-1]-wo,log(abs(fDIp(wo+pp,:))),'mo-','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/TIRFbasic/WoFilterCenter.m b/TIRFbasic/WoFilterCenter.m
new file mode 100644
index 0000000..8299473
--- /dev/null
+++ b/TIRFbasic/WoFilterCenter.m
@@ -0,0 +1,57 @@
+function [FiSMaof,NoisePower] = WoFilterCenter(FiSMao,OTFo,co,OBJpara,SFo)
+% Aim: Wiener Filtering the central frequency component
+% INPUT VARIABLES
+% FiSMao: noisy central frequency component
+% OTFo: system OTF
+% co: Wiener filter constant [=1, for minimum RMS estimate]
+% SFo: scaling factor (not significant here, so set to 1)
+% OUTPUT VARIABLES
+% FiSMaof: Wiener Filtered estimate of FiSMao
+% NoisePower: avg. noise power in FiSMao
+
+w = size(FiSMao,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Ro = sqrt( (X-wo).^2 + (Y-wo).^2 );
+
+OTFpower = OTFo.*conj(OTFo);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = FiSMao.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Object Power determination
+Aobj = OBJpara(1);
+Bobj = OBJpara(2);
+Ro(wo+1,wo+1) = 1;
+OBJpower = Aobj*(Ro.^Bobj);
+OBJpower = OBJpower.^2;
+
+%% Wiener Filtering
+FiSMaof = FiSMao.*(SFo.*conj(OTFo)./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+
+%% for cross-checking filtered estimate visually
+%{
+WFilter = (SFo.^2.*OTFpower./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+FiSMao1 = FiSMaof.*OTFo;
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(FiSMaof(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao(wo+pp,:))),'o--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao1(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,10*WFilter(wo+pp,:),'k*--','LineWidth',2,'MarkerSize',4)
+title('WoFilter')
+legend('FiSMaof','FiSMao','FiSMao1')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/TIRFbasic/WoFilterSideLobe.m b/TIRFbasic/WoFilterSideLobe.m
new file mode 100644
index 0000000..fb717c5
--- /dev/null
+++ b/TIRFbasic/WoFilterSideLobe.m
@@ -0,0 +1,54 @@
+function [FiSMaof,NoisePower] = WoFilterSideLobe(FiSMao,OTFo,co,OBJsideP,SFo)
+% Aim: Wiener Filtering the off-center frequency component
+% INPUT VARIABLES
+% FiSMao: noisy off-center frequency component
+% OTFo: system OTF
+% co: Wiener filter constant [=1, for minimum RMS estimate]
+% OBJsideP: avg. object spectrum
+% SFo: modulation factor
+% OUTPUT VARIABLES
+% FiSMaof: Wiener Filtered estimate of FiSMao
+% NoisePower: avg. noise power in FiSMao
+
+w = size(FiSMao,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Ro = sqrt( (X-wo).^2 + (Y-wo).^2 );
+
+OTFpower = OTFo.*conj(OTFo);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = FiSMao.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Object Power determination
+OBJpower = OBJsideP.^2;
+
+%% Wiener Filtering
+FiSMaof = FiSMao.*(SFo.*conj(OTFo)./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+
+%% for cross-checking filtered estimate visually
+%{
+WFilter = (SFo.^2.*OTFpower./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+FiSMao1 = FiSMaof.*OTFo;
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(FiSMaof(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao(wo+pp,:))),'o--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao1(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,10*WFilter(wo+pp,:),'k*--','LineWidth',2,'MarkerSize',4)
+title('WoFilter')
+legend('FiSMaof','FiSMao','FiSMao1')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/TIRFbasic/testpat.tiff b/TIRFbasic/testpat.tiff
new file mode 100644
index 0000000..84f1d62
Binary files /dev/null and b/TIRFbasic/testpat.tiff differ
diff --git a/TIRFexpt/ApodizationFunction.m b/TIRFexpt/ApodizationFunction.m
new file mode 100644
index 0000000..172f322
--- /dev/null
+++ b/TIRFexpt/ApodizationFunction.m
@@ -0,0 +1,30 @@
+function [FsumA] = ApodizationFunction(Fsum,k2fa,k2fb,k2fc,Kotf,Index)
+
+t = size(Fsum,1);
+
+OTFo_mask = OTFmaskShifted(0.*k2fa,Kotf,t);
+OTFap_mask = OTFmaskShifted(k2fa,Kotf,t);
+OTFam_mask = OTFmaskShifted(-k2fa,Kotf,t);
+OTFbp_mask = OTFmaskShifted(k2fb,Kotf,t);
+OTFbm_mask = OTFmaskShifted(-k2fb,Kotf,t);
+OTFcp_mask = OTFmaskShifted(k2fc,Kotf,t);
+OTFcm_mask = OTFmaskShifted(-k2fc,Kotf,t);
+ApoMask = OTFo_mask.*OTFap_mask.*OTFam_mask.*OTFbp_mask.*OTFbm_mask.*OTFcp_mask.*OTFcm_mask;
+clear OTFoz OTFo_mask OTFap_mask OTFam_mask OTFbp_mask OTFbm_mask OTFcp_mask OTFcm_mask
+%{
+figure;
+imshow(ApoMask,[ ])
+%}
+
+DistApoMask = bwdist(ApoMask);
+maxApoMask = max(max(DistApoMask));
+ApoFunc = double(DistApoMask./maxApoMask).^Index;
+%{
+figure;
+imshow(DistApoMask,[ ])
+figure;
+mesh(double(ApoFunc))
+%}
+
+FsumA = Fsum.*ApoFunc;
+clear ApoMask DistApoMask ApoFunc
\ No newline at end of file
diff --git a/TIRFexpt/BgNormalization0.m b/TIRFexpt/BgNormalization0.m
new file mode 100644
index 0000000..7d550b1
--- /dev/null
+++ b/TIRFexpt/BgNormalization0.m
@@ -0,0 +1,39 @@
+function [S1aTnoisyA, S2aTnoisyA, S3aTnoisyA, ...
+ S1bTnoisyA, S2bTnoisyA, S3bTnoisyA, ...
+ S1cTnoisyA, S2cTnoisyA, S3cTnoisyA] = BgNormalization0(...
+ S1aTnoisy, S2aTnoisy, S3aTnoisy, ...
+ S1bTnoisy, S2bTnoisy, S3bTnoisy, ...
+ S1cTnoisy, S2cTnoisy, S3cTnoisy)
+
+m1a = mean2(S1aTnoisy);
+m2a = mean2(S2aTnoisy);
+m3a = mean2(S3aTnoisy);
+m1b = mean2(S1bTnoisy);
+m2b = mean2(S2bTnoisy);
+m3b = mean2(S3bTnoisy);
+m1c = mean2(S1cTnoisy);
+m2c = mean2(S2cTnoisy);
+m3c = mean2(S3cTnoisy);
+
+s1a = std(S1aTnoisy(:));
+s2a = std(S2aTnoisy(:));
+s3a = std(S3aTnoisy(:));
+s1b = std(S1bTnoisy(:));
+s2b = std(S2bTnoisy(:));
+s3b = std(S3bTnoisy(:));
+s1c = std(S1cTnoisy(:));
+s2c = std(S2cTnoisy(:));
+s3c = std(S3cTnoisy(:));
+
+mTemp = max([m1a m2a m3a m1b m2b m3b m1c m2c m3c]);
+sTemp = max([s1a s2a s3a s1b s2b s3b s1c s2c s3c]);
+
+S1aTnoisyA = (S1aTnoisy - m1a).*(sTemp./s1a) + mTemp;
+S2aTnoisyA = (S2aTnoisy - m2a).*(sTemp./s2a) + mTemp;
+S3aTnoisyA = (S3aTnoisy - m3a).*(sTemp./s3a) + mTemp;
+S1bTnoisyA = (S1bTnoisy - m1b).*(sTemp./s1b) + mTemp;
+S2bTnoisyA = (S2bTnoisy - m2b).*(sTemp./s2b) + mTemp;
+S3bTnoisyA = (S3bTnoisy - m3b).*(sTemp./s3b) + mTemp;
+S1cTnoisyA = (S1cTnoisy - m1c).*(sTemp./s1c) + mTemp;
+S2cTnoisyA = (S2cTnoisy - m2c).*(sTemp./s2c) + mTemp;
+S3cTnoisyA = (S3cTnoisy - m3c).*(sTemp./s3c) + mTemp;
\ No newline at end of file
diff --git a/TIRFexpt/Ifreq2opt.m b/TIRFexpt/Ifreq2opt.m
new file mode 100644
index 0000000..02fe3bc
--- /dev/null
+++ b/TIRFexpt/Ifreq2opt.m
@@ -0,0 +1,56 @@
+function [CCop] = Ifreq2opt(k2fa,fAo,fAp,OTFo)
+% Aim: compute cross-correlation between central and off-central frequency
+% components
+% INPUT VARIABLES:
+% k2fa: nearest pixel approximation to illumination frequency vector
+% fAo: central frequency component
+% fAp: off-center frequency component
+% OTFo: system OTF
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% suppressing residual low frequency peaks of off-center components using
+% G is a notch-filter (determined heuristically)
+G = 1 - exp(-0.05*Ro.^1.2);
+fAo = fAo.*G;
+fAp = fAp.*G;
+
+fAoT = fAo.*conj(OTFo);
+fApT = fAp.*conj(OTFo);
+
+Kotf = OTFedgeF(OTFo);
+DoubleMatSize = 0;
+if ( 2*Kotf > wo )
+ DoubleMatSize = 1; % 1 for doubling fourier domain size, 0 for keeping it unchanged
+end
+
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ fAoT_temp = zeros(t,t);
+ fApT_temp = zeros(t,t);
+ fAoT_temp(wo+1:w+wo,wo+1:w+wo) = fAoT;
+ fApT_temp(wo+1:w+wo,wo+1:w+wo) = fApT;
+ clear fAoT fApT
+ fAoT = fAoT_temp;
+ fApT = fApT_temp;
+ clear fAoT_temp fApT_temp
+else
+ t = w;
+end
+to = t/2;
+u = linspace(0,t-1,t);
+v = linspace(0,t-1,t);
+[U,V] = meshgrid(u,v);
+ApT = exp( +1i.*2*pi*( k2fa(2)/t.*(U-to)+k2fa(1)/t.*(V-to) ) ).*ifft2(fApT);
+fApT0 = fft2( ApT );
+
+mA = sum(sum( fAoT.*conj(fApT0) ));
+mA = mA./sum(sum( fApT0.*conj(fApT0) ));
+
+CCop = -abs(mA);
\ No newline at end of file
diff --git a/TIRFexpt/IlluminationFreqTIRF.m b/TIRFexpt/IlluminationFreqTIRF.m
new file mode 100644
index 0000000..ec515c7
--- /dev/null
+++ b/TIRFexpt/IlluminationFreqTIRF.m
@@ -0,0 +1,71 @@
+function k2fa1 = IlluminationFreqTIRF(fDo,fDp,OTFo,k2o,thetaA)
+
+% AIM: To estimate illumination frequency vector
+% INPUT VARIABLES:
+% fDo: noisy central frequency component
+% fDp: noisy off-center frequency component
+% OTFo: system OTF
+% k2o: tentative magnitude of illumination frequency vector
+% thetaA: angular orientation of structured illumination
+
+%% for noise suppression
+fAp0 = fDp.*OTFo;
+fAo0 = fDo.*OTFo;
+
+% OTF cut-off freq
+Kotf = OTFedgeF(OTFo);
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+Zo = double( Ro<Kotf );
+
+Ka = k2o.*[cos(thetaA) sin(thetaA)];
+Ck = Ka(1) + 1i*Ka(2);
+rp = abs(Cv - Ck);
+rm = abs(Cv + Ck);
+r1 = 20; % local search radius
+Zmask = (rp>r1).*(rm>r1);
+Zmask = 1 - Zmask;
+
+% computing cross-correlation between central and off-central frequency
+% components to find nearest pixel approximation to illumination frequency
+% vector
+Kmax = k2o + r1;
+RHO = zeros(w,w);
+for p = wo-Kmax:wo+Kmax
+ for q = wo-Kmax:wo+Kmax
+ if ( Zmask(p,q) > 0 )
+ Cp = (p-wo) + 1i*(q-wo);
+ Rp = abs(Cv-Cp);
+ Zp = double( Rp<Kotf );
+ Z1 = Zo.*Zp;
+ fAp1 = circshift(fAp0,[(p-wo) (q-wo)]);
+ CC = sum(sum( fAo0.*conj(fAp1) ));
+ RHO(p,q) = CC./sum(sum(Z1));
+ end
+ end
+end
+% figure, mesh(abs(RHO))
+
+[Mx Idx] = max(abs(RHO),[],1);
+[Mx0 Idx0] = max(Mx);
+[My Idy] = max(abs(RHO),[],2);
+[My0 Idy0] = max(My);
+px0 = Idx0 - wo;
+py0 = Idy0 - wo;
+k2fa = [py0 px0]; % nearest pixel approximation
+
+%% subpixel approximation
+Ifreq2opt0 = @(k2fa0)Ifreq2opt(k2fa0,fDo,fDp,OTFo);
+% options = optimset('LargeScale','off','Algorithm',...
+% 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxIter',200,'Display','iter');
+k2fa0 = k2fa;
+[k2fa1,fval] = fminsearch(Ifreq2opt0,k2fa0,options);
+k2a = sqrt(k2fa1*k2fa1')
diff --git a/TIRFexpt/Kai2Opt.m b/TIRFexpt/Kai2Opt.m
new file mode 100644
index 0000000..9dae25a
--- /dev/null
+++ b/TIRFexpt/Kai2Opt.m
@@ -0,0 +1,37 @@
+function CCo = Kai2Opt(phaseShift,FcS1aT,FcS2aT,OTFo)
+% Aim: Computes normalized cross-power spectrum between separated frequency
+% components
+% INPUT VARIABLES:
+% phaseShift: relative phases [phaseShift2,phaseShift3]
+% FcS1aT: FT of difference between raw images 1 and 2
+% FcS2aT: FT of difference between raw images 1 and 3
+% OTFo: system OTF
+
+phaseShift2 = phaseShift(1);
+phaseShift3 = phaseShift(2);
+
+phaseShift1 = 0;
+phase2 = exp(-1i*phaseShift2) - exp(-1i*phaseShift1);
+phase3 = exp(-1i*phaseShift3) - exp(-1i*phaseShift1);
+%% Transformation Matrix
+M = [phase2 conj(phase2);
+ phase3 conj(phase3)];
+
+%% Separting the components
+%===========================================================
+Minv = inv(M);
+
+FiSMap = Minv(1,1)*FcS1aT + Minv(1,2)*FcS2aT;
+FiSMam = Minv(2,1)*FcS1aT + Minv(2,2)*FcS2aT;
+
+%{
+figure;
+mesh(log(abs(FiSMap)))
+figure;
+mesh(log(abs(FiSMam)))
+kk
+%}
+
+FiSMap0 = FiSMap.*conj(OTFo);
+FiSMam0 = FiSMam.*conj(OTFo);
+CCo = abs( sum(sum( FiSMap.*conj(FiSMam).*(1-0*OTFo.^10) )) );
diff --git a/TIRFexpt/MergingHeptaletsF.m b/TIRFexpt/MergingHeptaletsF.m
new file mode 100644
index 0000000..e62467d
--- /dev/null
+++ b/TIRFexpt/MergingHeptaletsF.m
@@ -0,0 +1,65 @@
+function [Fsum,Fperi,Fcent] = MergingHeptaletsF(fDIo,fDIp,fDIm,fDBo,fDBp,fDBm,...
+ fDCo,fDCp,fDCm,Ma,Mb,Mc,npDIo,npDIp,npDIm,...
+ npDBo,npDBp,npDBm,npDCo,npDCp,npDCm,k2fa,k2fb,k2fc,OBJpara,OTFo)
+
+% AIM: To merge all 9 frequency components into one using generalized
+% Wiener Filter
+% INPUT VARIABLES
+% [fDIo fDIp fDIm
+% fDBo fDBp fDBm
+% fDCo fDCp fDCm]: nine frequency components
+% Ma,Mb,Mc: modulation factors for the three illumination orientations
+% [npDIo npDIp npDIm
+% npDBo npDBp npDBm
+% npDCo,npDCp,npDCm]: noise powers corresponding to nine frequency
+% components
+% k2fa,k2fb,k2fc: illumination frequency vectors for the three
+% illumination orientations
+% OBJpara: Object spectrum parameters
+% OTFo: system OTF
+% OUTPUT VARIABLES
+% Fsum: all nine frequency components merged into one using
+% generalised Wiener Filter
+% Fperi: six off-center frequency components merged into one using
+% generalised Wiener Filter
+% Fcent: averaged of the three central frequency components
+
+% obtain signal spectrums corresponding to central and off-center
+% frequency components
+[SIGao,SIGam2,SIGap2] = TripletSNR0(OBJpara,k2fa,OTFo,fDIm);
+[SIGbo,SIGbm2,SIGbp2] = TripletSNR0(OBJpara,k2fb,OTFo,fDBm);
+[SIGco,SIGcm2,SIGcp2] = TripletSNR0(OBJpara,k2fc,OTFo,fDCm);
+
+SIGap2 = Ma*SIGap2;
+SIGam2 = Ma*SIGam2;
+SIGbp2 = Mb*SIGbp2;
+SIGbm2 = Mb*SIGbm2;
+SIGcp2 = Mc*SIGcp2;
+SIGcm2 = Mc*SIGcm2;
+
+%% Generalized Wiener-Filter computation
+SNRao = SIGao.*conj(SIGao)./npDIo;
+SNRap = SIGap2.*conj(SIGap2)./npDIp;
+SNRam = SIGam2.*conj(SIGam2)./npDIm;
+
+SNRbo = SIGbo.*conj(SIGbo)./npDBo;
+SNRbp = SIGbp2.*conj(SIGbp2)./npDBp;
+SNRbm = SIGbm2.*conj(SIGbm2)./npDBm;
+
+SNRco = SIGco.*conj(SIGco)./npDCo;
+SNRcp = SIGcp2.*conj(SIGcp2)./npDCp;
+SNRcm = SIGcm2.*conj(SIGcm2)./npDCm;
+
+ComDeno = 0.01 + ( SNRao + SNRap + SNRam + SNRao + SNRbp + SNRbm + SNRco + SNRcp + SNRcm );
+Fsum = fDIo.*SNRao + fDIp.*SNRap + fDIm.*SNRam...
+ + fDBo.*SNRbo + fDBp.*SNRbp + fDBm.*SNRbm...
+ + fDCo.*SNRco + fDCp.*SNRcp + fDCm.*SNRcm;
+Fsum = Fsum./ComDeno;
+
+ComPeri = 0.01 + ( SNRap + SNRam + SNRbp + SNRbm + SNRcp + SNRcm );
+Fperi = fDIp.*SNRap + fDIm.*SNRam + fDBp.*SNRbp + fDBm.*SNRbm...
+ + fDCp.*SNRcp + fDCm.*SNRcm;
+Fperi = Fperi./ComPeri;
+
+% averaged central frequency component
+Fcent = (fDIo + fDBo + fDCo)/3;
diff --git a/TIRFexpt/ModulationFactor.m b/TIRFexpt/ModulationFactor.m
new file mode 100644
index 0000000..2ac2dda
--- /dev/null
+++ b/TIRFexpt/ModulationFactor.m
@@ -0,0 +1,86 @@
+function Mm = ModulationFactor(fDp,k1a,OBJparaA,OTFo)
+% AIM: Determination of modulation factor
+% INPUT VARIABLES
+% fDp: off-center frequency component
+% k1a: illumination frequency vector
+% OBJparaA: Object power parameters
+% OTFo: system OTF
+% OUTPUT VARIABLE
+% Mm: modulation factor
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% magnitude of illumination vector
+k2 = sqrt(k1a*k1a');
+
+% vector along illumination direction
+kv = k1a(2) + 1i*k1a(1);
+Rp = abs(Cv+kv);
+
+% Object power parameters
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+
+% Object spectrum
+OBJp = Aobj*(Rp+0).^Bobj;
+
+% illumination vector rounded to nearest pixel
+k3 = -round(k1a);
+
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+
+% signal spectrum
+SIGap = OBJp.*OTFo;
+
+% figure;
+% mesh(log(abs(OBJp)))
+% figure;
+% mesh(log(abs(fDp)))
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = fDp.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Noise free object power computation
+Fpower = fDp.*conj(fDp) - NoisePower;
+fDp = sqrt(abs(Fpower));
+
+% frequency range over which signal power matching is done to estimate
+% modulation factor
+Zmask = (Ro > 0.2*k2).*(Ro < 0.8*k2).*(Rp > 0.2*k2);
+
+% least square approximation for modulation factor
+Mm = sum(sum(SIGap.*abs(fDp).*Zmask));
+Mm = Mm./sum(sum(SIGap.^2.*Zmask))
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(fDp(wo+pp,:))),'o-','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(SIGap(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,log(abs(Mm*SIGap(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+grid on
+box on
+kk
+%}
+
+
+
diff --git a/TIRFexpt/OBJparaOpt.m b/TIRFexpt/OBJparaOpt.m
new file mode 100644
index 0000000..8acbbbb
--- /dev/null
+++ b/TIRFexpt/OBJparaOpt.m
@@ -0,0 +1,47 @@
+function Esum = OBJparaOpt(OBJpara0,fDIoTnoisy,OTFo)
+% AIM: Determined Sum of Squared Errors (SSE) between `actual signal power' and
+% `approximated signal power'
+% INPUT VARIABLES
+% OBJpara0: [Aobj Bobj], Object power parameters
+% fDIoTnoisy: FT of central frequency component
+% OTFo: system OTF
+% OUTPUT VARIABLE
+% Esum: SSE between `actual signal spectrum' and `approximated signal spectrum'
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+Ro(wo+1,wo+1) = 1; % to avoid nan
+
+% approximated signal power calculation
+Aobj = OBJpara0(1);
+Bobj = OBJpara0(2);
+OBJpower = Aobj*(Ro.^Bobj);
+SIGpower = OBJpower.*OTFo;
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% range of frequency over which SSE is computed
+Zloop = (Ro<0.75*Kotf).*(Ro>0.25*Kotf);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = fDIoTnoisy.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Noise free object power computation
+Fpower = fDIoTnoisy.*conj(fDIoTnoisy) - NoisePower;
+fDIoTnoisy = sqrt(abs(Fpower));
+
+% SSE computation
+Error = fDIoTnoisy - SIGpower;
+Esum = sum(sum((Error.^2./Ro).*Zloop));
+
diff --git a/TIRFexpt/OBJpowerPara.m b/TIRFexpt/OBJpowerPara.m
new file mode 100644
index 0000000..9ef86ce
--- /dev/null
+++ b/TIRFexpt/OBJpowerPara.m
@@ -0,0 +1,48 @@
+function OBJparaA = OBJpowerPara(fDIoTnoisy,OTFo)
+% AIM: determination of object power parameters Aobj and Bobj
+% INPUT VARIABLES
+% fDIoTnoisy: FT of central frequency component
+% OTFo: system OTFo
+% OUTPUT VARIABLES
+% OBJparaA: [Aobj Bobj]
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+%% object power parameters through optimization
+OBJparaOpt0 = @(OBJpara0)OBJparaOpt(OBJpara0,fDIoTnoisy,OTFo);
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+
+% obtaining crude initial guesses for Aobj and Bobj
+Zm = (Ro>0.3*Kotf).*(Ro<0.4*Kotf);
+Aobj = sum(sum(abs(fDIoTnoisy.*Zm)))./sum(sum(Zm));
+Bobj = -0.5;
+OBJpara0 = [Aobj Bobj];
+
+% optimization step
+[OBJparaA,fval] = fminsearch(OBJparaOpt0,OBJpara0,options);
+
+%% plot for cross-checking the result
+%{
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+OBJo = Aobj*(Ro.^Bobj);
+SIGp = OBJo.*OTFo;
+pp = 3;
+figure;
+hold on
+plot([0:w-1]-wo,log(abs(fDIoTnoisy(wo+pp,:))),'k--','LineWidth',3,'MarkerSize',6)
+plot([0:w-1]-wo,log(abs(SIGp(wo+pp,:))),'mo-','LineWidth',2,'MarkerSize',6)
+legend('acutal signal power','avg. signal power')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/TIRFexpt/OTF.tif b/TIRFexpt/OTF.tif
new file mode 100644
index 0000000..424998d
Binary files /dev/null and b/TIRFexpt/OTF.tif differ
diff --git a/TIRFexpt/OTFdoubling.m b/TIRFexpt/OTFdoubling.m
new file mode 100644
index 0000000..58aaa6e
--- /dev/null
+++ b/TIRFexpt/OTFdoubling.m
@@ -0,0 +1,17 @@
+function OTF2 = OTFdoubling(OTFo,DoubleMatSize)
+
+%% embeds OTFo in doubled range of frequencies
+
+w = size(OTFo,1);
+wo = w/2;
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ OTF2 = zeros(2*w,2*w);
+ OTF2(wo+1:w+wo,wo+1:w+wo) = OTFo;
+ clear OTFo
+else
+ OTF2 = OTFo;
+ clear OTFo
+end
\ No newline at end of file
diff --git a/TIRFexpt/OTFedgeF.m b/TIRFexpt/OTFedgeF.m
new file mode 100644
index 0000000..e4d6a7f
--- /dev/null
+++ b/TIRFexpt/OTFedgeF.m
@@ -0,0 +1,13 @@
+function Kotf = OTFedgeF(OTFo)
+
+w = size(OTFo,1);
+wo = w/2;
+
+OTF1 = OTFo(wo+1,:);
+OTFmax = max(max(abs(OTFo)));
+OTFtruncate = 0.01;
+i = 1;
+while ( abs(OTF1(1,i))<OTFtruncate*OTFmax )
+ Kotf = wo+1-i;
+ i = i + 1;
+end
\ No newline at end of file
diff --git a/TIRFexpt/OTFmaskShifted.m b/TIRFexpt/OTFmaskShifted.m
new file mode 100644
index 0000000..a60ac98
--- /dev/null
+++ b/TIRFexpt/OTFmaskShifted.m
@@ -0,0 +1,13 @@
+function [OTFap_mask] = OTFmaskShifted(k2fa,Kotf,w)
+% OTFap_mask = complement of shifted OTF
+
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+
+Iy = k2fa(1,1) + (wo+1);
+Ix = k2fa(1,2) + (wo+1);
+
+Ro = sqrt( (X+1-Ix).^2 + (Y+1-Iy).^2 );
+OTFap_mask = Ro > Kotf;
\ No newline at end of file
diff --git a/TIRFexpt/OTFpost.m b/TIRFexpt/OTFpost.m
new file mode 100644
index 0000000..2498868
--- /dev/null
+++ b/TIRFexpt/OTFpost.m
@@ -0,0 +1,20 @@
+function OTFo = OTFpost(OTFo)
+%% AIM: To flaten out the peripheral rippling frequencies
+
+OTFo = OTFo./max(OTFo(:));
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+Zmask = (Ro>(Kotf+20)).*(Ro<wo);
+OTFperi = sum(sum(OTFo.*Zmask))./sum(sum(Zmask));
+Zperi = (Ro>(wo-10));
+OTFo = OTFo.*(1-Zperi) + OTFperi.*Zperi;
\ No newline at end of file
diff --git a/TIRFexpt/PCMfilteringF.m b/TIRFexpt/PCMfilteringF.m
new file mode 100644
index 0000000..19d944a
--- /dev/null
+++ b/TIRFexpt/PCMfilteringF.m
@@ -0,0 +1,172 @@
+function [fDof,fDp2,fDm2,npDo,npDp,npDm,Mm,DoubleMatSize]...
+ = PCMfilteringF(fDo,fDp,fDm,OTFo,OBJparaA,kA)
+
+% AIM: obtaining Wiener Filtered estimates of noisy frequency components
+% INPUT VARIABLES
+% fDo,fDp,fDm: noisy estimates of separated frequency components
+% OTFo: system OTF
+% OBJparaA: object power parameters
+% kA: illumination vector
+% OUTPUT VARIABLES
+% fDof,fDp2,fDm2: Wiener Filtered estimates of fDo,fDp,fDm
+% npDo,npDp,npDm: avg. noise power in fDo,fDp,fDm
+% Mm: illumination modulation factor
+% DoubleMatSize: parameter for doubling FT size if necessary
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% suppressing out-of-focus signal of off-center components using
+% notch-filter (determined heuristically)
+G = 1 - exp(-0.05*Ro.^1.2);
+fDp = fDp.*G;
+fDm = fDm.*G;
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% object power parameters
+Aobj = OBJparaA(1);
+Bobj = OBJparaA(2);
+
+% Wiener Filtering central frequency component
+SFo = 1;
+co = 1.0;
+[fDof,npDo] = WoFilterCenter(fDo,OTFo,co,OBJparaA,SFo);
+
+% modulation factor determination
+Mm = ModulationFactor(fDp,kA,OBJparaA,OTFo);
+
+%% Duplex power (default)
+kv = kA(2) + 1i*kA(1); % vector along illumination direction
+Rp = abs(Cv-kv);
+Rm = abs(Cv+kv);
+OBJp = Aobj*(Rp.^Bobj);
+OBJm = Aobj*(Rm.^Bobj);
+k3 = round(kA);
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+OBJm(wo+1-k3(1),wo+1-k3(2)) = 0.25*OBJm(wo+2-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+2-k3(2))...
+ + 0.25*OBJm(wo+0-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+0-k3(2));
+
+% Filtering side lobes (off-center frequency components)
+SFo = Mm;
+[fDpf,npDp] = WoFilterSideLobe(fDp,OTFo,co,OBJm,SFo);
+[fDmf,npDm] = WoFilterSideLobe(fDm,OTFo,co,OBJp,SFo);
+
+%% doubling Fourier domain size if necessary
+DoubleMatSize = 0;
+if ( 2*Kotf > wo )
+ DoubleMatSize = 1; % 1 for doubling fourier domain size, 0 for keeping it unchanged
+end
+if ( DoubleMatSize>0 )
+ t = 2*w;
+ to = t/2;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ [U,V] = meshgrid(u,v);
+ fDoTemp = zeros(2*w,2*w);
+ fDpTemp = zeros(2*w,2*w);
+ fDmTemp = zeros(2*w,2*w);
+ OTFtemp = zeros(2*w,2*w);
+ fDoTemp(wo+1:w+wo,wo+1:w+wo) = fDof;
+ fDpTemp(wo+1:w+wo,wo+1:w+wo) = fDpf;
+ fDmTemp(wo+1:w+wo,wo+1:w+wo) = fDmf;
+ OTFtemp(wo+1:w+wo,wo+1:w+wo) = OTFo;
+ clear fDof fDpf fDmf OTFo w wo x y X Y
+ fDof = fDoTemp;
+ fDpf = fDpTemp;
+ fDmf = fDmTemp;
+ OTFo = OTFtemp;
+ clear fDoTemp fDpTemp fDmTemp OTFtemp
+else
+ t = w;
+ to = t/2;
+ u = linspace(0,t-1,t);
+ v = linspace(0,t-1,t);
+ [U,V] = meshgrid(u,v);
+end
+
+% Shifting the off-center frequency components to their correct location
+fDp1 = fft2(ifft2(fDpf).*exp( +1i.*2*pi*(kA(2)/t.*(U-to) + kA(1)/t.*(V-to)) ));
+fDm1 = fft2(ifft2(fDmf).*exp( -1i.*2*pi*(kA(2)/t.*(U-to) + kA(1)/t.*(V-to)) ));
+%{
+figure;
+mesh(log(abs(fDpf)))
+figure;
+mesh(log(abs(fDp1)))
+kk
+%}
+
+%% Shift induced phase error correction
+Cv = (U-to) + 1i*(V-to);
+Ro = abs(Cv);
+Rp = abs(Cv-kv);
+k2 = sqrt(kA*kA');
+
+% frequency range over which corrective phase is determined
+Zmask = (Ro < 0.8*k2).*(Rp < 0.8*k2);
+
+% corrective phase
+Angle0 = angle( sum(sum( fDof.*conj(fDp1).*Zmask )) );
+
+% phase correction
+fDp2 = exp(+1i*Angle0).*fDp1;
+fDm2 = exp(-1i*Angle0).*fDm1;
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot(u-to,angle(fDof(to+pp,:)).*180/pi,'o-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(fDp1(to+pp,:)).*180/pi,'ro-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(fDp2(to+pp,:)).*180/pi,'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+figure;
+hold on
+plot3(u-to,real(fDof(to+1,:)),imag(fDof(to+1,:)),'o-','LineWidth',2,'MarkerSize',8)
+plot3(u-to,real(fDp1(to+1,:)),imag(fDp1(to+1,:)),'ro-','LineWidth',2,'MarkerSize',8)
+plot3(u-to,real(fDp2(to+1,:)),imag(fDp2(to+1,:)),'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+
+ff1 = fDof.*conj(fDp1);
+ff2 = fDof.*conj(fDp2);
+figure;
+hold on
+plot(u-to,angle(ff1(to+pp,:)).*180/pi,'ro-','LineWidth',2,'MarkerSize',8)
+plot(u-to,angle(ff2(to+pp,:)).*180/pi,'go--','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+%}
+
+%{
+fperi = fDp2 + fDm2;
+fall = fDof + fperi;
+
+Dof = real( ifft2(fftshift(fDof)) );
+figure;
+imshow(Dof,[])
+
+Dperi = real( ifft2(fftshift(fperi)) );
+figure;
+imshow(Dperi,[])
+
+Dall = real( ifft2(fftshift(fall)) );
+figure;
+imshow(Dall,[])
+%}
+
+
+
diff --git a/TIRFexpt/PCMseparateF.m b/TIRFexpt/PCMseparateF.m
new file mode 100644
index 0000000..70e2d56
--- /dev/null
+++ b/TIRFexpt/PCMseparateF.m
@@ -0,0 +1,57 @@
+function [fDo,fDp,fDm]...
+ = PCMseparateF(S1aTnoisy,S2aTnoisy,S3aTnoisy,OTFo)
+% AIM: obtaining the noisy estimates of three frequency components
+% INPUT VARIABLES
+% S1aTnoisy,S2aTnoisy,S3aTnoisy: 3 raw SIM images with identical
+% illumination pattern orientation
+% but different phase shifts
+% OTFo: system OTF
+% OUTPUT VARIABLES
+% fDo,fDp,fDm: noisy estimates of separated frequency components
+
+
+w = size(S1aTnoisy,1);
+wo = w/2;
+
+% computing PSFe for edge tapering SIM images
+PSFd = real(fftshift( ifft2(fftshift(OTFo.^3)) ));
+PSFd = PSFd/max(max(PSFd));
+PSFd = PSFd/sum(sum(PSFd));
+h = 30;
+PSFe = PSFd(wo-h+1:wo+h,wo-h+1:wo+h);
+
+% edge tapering raw SIM images
+S1aTnoisy = edgetaper(S1aTnoisy,PSFe);
+S2aTnoisy = edgetaper(S2aTnoisy,PSFe);
+S3aTnoisy = edgetaper(S3aTnoisy,PSFe);
+
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+fS2aTnoisy = fftshift(fft2(S2aTnoisy));
+fS3aTnoisy = fftshift(fft2(S3aTnoisy));
+
+fDuplex2 = fS2aTnoisy - fS1aTnoisy;
+fDuplex3 = fS3aTnoisy - fS1aTnoisy;
+
+%% Optimizing the relative phases
+Kai2Opt0 = @(phase0)Kai2Opt(phase0,fDuplex2,fDuplex3,OTFo);
+options = optimset('LargeScale','off','Algorithm',...
+ 'active-set','MaxFunEvals',500,'MaxIter',500,'Display','notify');
+phase0 = [pi/2 -pi/2]; % initial guess
+[phaseShift,fval] = fminsearch(Kai2Opt0,phase0,options);
+phaseShift*180/pi
+
+%% Separating the three frequency components
+phaseShift0 = 0;
+[fDo,fDp,fDm] = SeparatedComponents2D(...
+ phaseShift,phaseShift0,fS1aTnoisy,fS2aTnoisy,fS3aTnoisy);
+
+%% for visual verification
+%{
+figure;
+mesh(log(abs(fDo)))
+figure;
+mesh(log(abs(fDp)))
+figure;
+mesh(log(abs(fDm)))
+% kk
+%}
\ No newline at end of file
diff --git a/TIRFexpt/SIMo.m b/TIRFexpt/SIMo.m
new file mode 100644
index 0000000..d69e1ad
--- /dev/null
+++ b/TIRFexpt/SIMo.m
@@ -0,0 +1,94 @@
+clear all
+close all
+clc
+
+%% Loading system OTF file
+OTFo = double(imread('OTF.tif'));
+OTFo = OTFpost(OTFo);
+
+%% Read Expt. Raw SIM Images
+aa1 = imread('sim01z4.tif');
+aa = aa1(:,:,2);
+bb = uint16(zeros(512,512,9));
+for ii=1:3
+ for jj=1:3
+ bb(1:512,1:512,(ii-1)*3+jj)=aa((ii-1)*512+1:ii*512,(jj-1)*512+1:jj*512,1);
+ end
+end
+
+% separating the raw SIM images
+S1aTnoisy = double( bb(:,:,1) );
+S2aTnoisy = double( bb(:,:,2) );
+S3aTnoisy = double( bb(:,:,3) );
+S1bTnoisy = double( bb(:,:,4) );
+S2bTnoisy = double( bb(:,:,5) );
+S3bTnoisy = double( bb(:,:,6) );
+S1cTnoisy = double( bb(:,:,7) );
+S2cTnoisy = double( bb(:,:,8) );
+S3cTnoisy = double( bb(:,:,9) );
+clear aa aa1 bb
+
+%% Intensity Normalization
+[S1aTnoisy, S2aTnoisy, S3aTnoisy, ...
+ S1bTnoisy, S2bTnoisy, S3bTnoisy, ...
+ S1cTnoisy, S2cTnoisy, S3cTnoisy] = BgNormalization0(...
+ S1aTnoisy, S2aTnoisy, S3aTnoisy, ...
+ S1bTnoisy, S2bTnoisy, S3bTnoisy, ...
+ S1cTnoisy, S2cTnoisy, S3cTnoisy);
+
+%% Background subtraction
+SE = strel('disk',20);
+S1aTnoisy = S1aTnoisy - imopen(S1aTnoisy,SE);
+S2aTnoisy = S2aTnoisy - imopen(S2aTnoisy,SE);
+S3aTnoisy = S3aTnoisy - imopen(S3aTnoisy,SE);
+S1bTnoisy = S1bTnoisy - imopen(S1bTnoisy,SE);
+S2bTnoisy = S2bTnoisy - imopen(S2bTnoisy,SE);
+S3bTnoisy = S3bTnoisy - imopen(S3bTnoisy,SE);
+S1cTnoisy = S1cTnoisy - imopen(S1cTnoisy,SE);
+S2cTnoisy = S2cTnoisy - imopen(S2cTnoisy,SE);
+S3cTnoisy = S3cTnoisy - imopen(S3cTnoisy,SE);
+
+% tentative values to define search region for illumination frequency
+% determination
+k2o = 125; % tentative illumination frequency
+thetaA = 20*pi/180; % orientations of structured illumination
+thetaB = 140*pi/180;
+thetaC = 260*pi/180;
+
+%% obtaining the noisy estimates of three frequency components
+[fAo,fAp,fAm]...
+ = PCMseparateF(S1aTnoisy,S2aTnoisy,S3aTnoisy,OTFo);
+kA = IlluminationFreqTIRF(fAo,fAp,OTFo,k2o,thetaA);
+[fBo,fBp,fBm]...
+ = PCMseparateF(S1bTnoisy,S2bTnoisy,S3bTnoisy,OTFo);
+kB = IlluminationFreqTIRF(fBo,fBp,OTFo,k2o,thetaB);
+[fCo,fCp,fCm]...
+ = PCMseparateF(S1cTnoisy,S2cTnoisy,S3cTnoisy,OTFo);
+kC = IlluminationFreqTIRF(fCo,fCp,OTFo,k2o,thetaC);
+
+% averaging the central frequency components
+fCent = (fAo + fBo + fCo)/3;
+
+% Object power parameters determination
+OBJparaA = OBJpowerPara(fCent,OTFo);
+
+%% Wiener Filtering the noisy frequency components
+[fAof,fApf,fAmf,Nao,Nap,Nam,Ma,DoubleMatSize]...
+ = PCMfilteringF(fAo,fAp,fAm,OTFo,OBJparaA,kA);
+[fBof,fBpf,fBmf,Nbo,Nbp,Nbm,Mb,~]...
+ = PCMfilteringF(fAo,fBp,fBm,OTFo,OBJparaA,kB);
+[fCof,fCpf,fCmf,Nco,Ncp,Ncm,Mc,~]...
+ = PCMfilteringF(fAo,fCp,fCm,OTFo,OBJparaA,kC);
+
+%% doubling Fourier domain size if necessary
+OTFo = OTFdoubling(OTFo,DoubleMatSize);
+
+%% merging all 9 frequency components using generalized Wiener Filter
+[Fsum,Fperi,Fcent] = MergingHeptaletsF(fAof,fApf,fAmf,...
+ fBof,fBpf,fBmf,fCof,fCpf,fCmf,...
+ Ma,Mb,Mc,Nao,Nap,Nam,Nbo,Nbp,Nbm,...
+ Nco,Ncp,Ncm,kA,kB,kC,OBJparaA,OTFo);
+
+% Plotting SIM results
+SIMplot(Fsum,Fperi,Fcent,OTFo,kA,kB,kC,S1aTnoisy);
+
diff --git a/TIRFexpt/SIMplot.m b/TIRFexpt/SIMplot.m
new file mode 100644
index 0000000..b83fd74
--- /dev/null
+++ b/TIRFexpt/SIMplot.m
@@ -0,0 +1,84 @@
+function [] = SIMplot(Fsum,Fperi,Fcent,OTFo,kA,kB,kC,S1aTnoisy)
+
+%% recontructed SIM images
+Dcent = real( ifft2(fftshift(Fcent)) );
+Dsum = real( ifft2(fftshift(Fsum)) );
+Dperi = real( ifft2(fftshift(Fperi)) );
+
+t = size(OTFo,1);
+h = 1*30; % for removing the image edges
+figure;
+imshow(Dcent(h+1:t-h,h+1:t-h),[])
+title('Wiener-Filtered wide-field')
+figure;
+imshow(Dsum(h+1:t-h,h+1:t-h),[])
+title('SIM image')
+figure;
+imshow(Dperi(h+1:t-h,h+1:t-h),[])
+title('SIM image (using only off-center frequency components)')
+
+% appodizing the merged frequency components
+Index = 0.4;
+Kotf = OTFedgeF(OTFo);
+[FsumA] = ApodizationFunction(Fsum,kA,kB,kC,Kotf,Index);
+[FperiA] = ApodizationFunction(Fperi,kA,kB,kC,Kotf,Index);
+DsumA = real( ifft2(fftshift(FsumA)) );
+DperiA = real( ifft2(fftshift(FperiA)) );
+
+figure;
+imshow(DsumA(h+1:t-h,h+1:t-h),[])
+title('appodized SIM image')
+figure;
+imshow(DperiA(h+1:t-h,h+1:t-h),[])
+title('appodized SIM image (using only off-center frequency components)')
+
+
+%% Frequency Plots
+fS1aTnoisy = fftshift(fft2(S1aTnoisy));
+w = size(OTFo,1);
+wo = w/2;
+w1 = size(fS1aTnoisy,1);
+if ( w > w1 )
+ fS1aTnoisy0 = zeros(w,w);
+ fS1aTnoisy0(wo-w1/2+1:wo+w1/2,wo-w1/2+1:wo+w1/2) = fS1aTnoisy;
+ clear fS1aTnoisy
+ fS1aTnoisy = fS1aTnoisy0;
+ clear fS1aTnoisy0;
+end
+
+
+p = 10;
+minL1 = min(min( abs(fS1aTnoisy).^(1/p) ));
+minL2 = min(min( abs(Fcent).^(1/p) ));
+minL3 = min(min( abs(Fsum).^(1/p) ));
+minL4 = min(min( abs(Fperi).^(1/p) ));
+minL5 = min(min( abs(FsumA).^(1/p) ));
+maxL1 = max(max( abs(fS1aTnoisy).^(1/p) ));
+maxL2 = max(max( abs(Fcent).^(1/p) ));
+maxL3 = max(max( abs(Fsum).^(1/p) ));
+maxL4 = max(max( abs(Fperi).^(1/p) ));
+maxL5 = max(max( abs(FsumA).^(1/p) ));
+minL = min([minL1,minL2,minL3,minL4,minL5]);
+maxL = max([maxL1,maxL2,maxL3,maxL4,maxL5]);
+
+figure;
+imshow(abs(fS1aTnoisy).^(1/p),[minL maxL])
+title('fS1aTnoisy')
+
+figure;
+imshow(abs(Fcent).^(1/p),[minL maxL])
+title('Weiner Filtered frequency')
+figure;
+imshow(abs(Fsum).^(1/p),[minL maxL])
+title('SIM frequency')
+figure;
+imshow(abs(Fperi).^(1/p),[minL maxL])
+title('SIM (off-center frequency components)')
+figure;
+imshow(abs(FsumA).^(1/p),[minL maxL])
+title('appodized SIM frequency')
+figure;
+imshow(abs(FperiA).^(1/p),[minL maxL])
+title('appodized SIM (off-center frequency components)')
+
+
diff --git a/TIRFexpt/SeparatedComponents2D.m b/TIRFexpt/SeparatedComponents2D.m
new file mode 100644
index 0000000..d124752
--- /dev/null
+++ b/TIRFexpt/SeparatedComponents2D.m
@@ -0,0 +1,22 @@
+function [FiSMao,FiSMap,FiSMam] = SeparatedComponents2D(...
+ phaseShift,phaseShift0,FcS1aT,FcS2aT,FcS3aT)
+% Aim: Unmixing the frequency components of raw SIM images
+% phaseShift,phaseShift0: illumination phase shifts
+% FcS1aT,FcS2aT,FcS3aT: FT of raw SIM images
+% FiSMao,FiSMap,FiSMam: unmixed frequency components of raw SIM images
+
+phaseShift1 = phaseShift(1);
+phaseShift2 = phaseShift(2);
+MF = 1.0;
+%% Transformation Matrix
+M = 0.5*[1 0.5*MF*exp(-1i*phaseShift0) 0.5*MF*exp(+1i*phaseShift0);
+ 1 0.5*MF*exp(-1i*phaseShift1) 0.5*MF*exp(+1i*phaseShift1);
+ 1 0.5*MF*exp(-1i*phaseShift2) 0.5*MF*exp(+1i*phaseShift2)];
+
+%% Separting the components
+%===========================================================
+Minv = inv(M);
+
+FiSMao = Minv(1,1)*FcS1aT + Minv(1,2)*FcS2aT + Minv(1,3)*FcS3aT;
+FiSMap = Minv(2,1)*FcS1aT + Minv(2,2)*FcS2aT + Minv(2,3)*FcS3aT;
+FiSMam = Minv(3,1)*FcS1aT + Minv(3,2)*FcS2aT + Minv(3,3)*FcS3aT;
diff --git a/TIRFexpt/TripletSNR0.m b/TIRFexpt/TripletSNR0.m
new file mode 100644
index 0000000..6bf5b1d
--- /dev/null
+++ b/TIRFexpt/TripletSNR0.m
@@ -0,0 +1,66 @@
+function [SIGao,SIGap2,SIGam2] = TripletSNR0(OBJpara,k2fa,OTFo,fDIp)
+
+% AIM: To obtain signal spectrums corresponding to central and off-center
+% frequency components
+% INPUT VARIABLES
+% OBJpara: object power parameters
+% k2fa: illumination frequency vector
+% OTFo: system OTF
+% fDIp: one of the off-center frequency component (utilized here only for
+% visual verification of computation)
+% OUTPUT VARIABLES
+% SIGao: signal spectrum corresponding to central frequency component
+% SIGap2,SIGam2: signal spectrums corresponding to off-center frequency components
+
+
+w = size(OTFo,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Cv = (X-wo) + 1i*(Y-wo);
+Ro = abs(Cv);
+
+% object spectrum parameters
+Aobj = OBJpara(1);
+Bobj = OBJpara(2);
+
+% object spectrums (default)
+kv = k2fa(2) + 1i*k2fa(1); % vector along illumination direction
+Rp = abs(Cv-kv);
+Rm = abs(Cv+kv);
+OBJo = Aobj*(Ro.^Bobj);
+OBJp = Aobj*(Rp.^Bobj);
+OBJm = Aobj*(Rm.^Bobj);
+
+OBJo(wo+1,wo+1) = 0.25*OBJo(wo+2,wo+1) + 0.25*OBJo(wo+1,wo+2)...
+ + 0.25*OBJo(wo+0,wo+1) + 0.25*OBJo(wo+1,wo+0);
+
+k3 = round(k2fa);
+OBJp(wo+1+k3(1),wo+1+k3(2)) = 0.25*OBJp(wo+2+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+2+k3(2))...
+ + 0.25*OBJp(wo+0+k3(1),wo+1+k3(2))...
+ + 0.25*OBJp(wo+1+k3(1),wo+0+k3(2));
+OBJm(wo+1-k3(1),wo+1-k3(2)) = 0.25*OBJm(wo+2-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+2-k3(2))...
+ + 0.25*OBJm(wo+0-k3(1),wo+1-k3(2))...
+ + 0.25*OBJm(wo+1-k3(1),wo+0-k3(2));
+
+% signal spectrums
+SIGao = OBJo.*OTFo;
+SIGap = OBJp.*OTFo;
+SIGam = OBJm.*OTFo;
+
+SIGap2 = circshift(SIGap,-k3);
+SIGam2 = circshift(SIGam,k3);
+
+%% for visual verification
+%{
+pp = 3;
+figure;
+hold on
+plot([0:w-1]-wo,log(abs(SIGap2(wo+pp,:))),'k--','LineWidth',3,'MarkerSize',6)
+plot([0:w-1]-wo,log(abs(fDIp(wo+pp,:))),'mo-','LineWidth',2,'MarkerSize',6)
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/TIRFexpt/WoFilterCenter.m b/TIRFexpt/WoFilterCenter.m
new file mode 100644
index 0000000..8299473
--- /dev/null
+++ b/TIRFexpt/WoFilterCenter.m
@@ -0,0 +1,57 @@
+function [FiSMaof,NoisePower] = WoFilterCenter(FiSMao,OTFo,co,OBJpara,SFo)
+% Aim: Wiener Filtering the central frequency component
+% INPUT VARIABLES
+% FiSMao: noisy central frequency component
+% OTFo: system OTF
+% co: Wiener filter constant [=1, for minimum RMS estimate]
+% SFo: scaling factor (not significant here, so set to 1)
+% OUTPUT VARIABLES
+% FiSMaof: Wiener Filtered estimate of FiSMao
+% NoisePower: avg. noise power in FiSMao
+
+w = size(FiSMao,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Ro = sqrt( (X-wo).^2 + (Y-wo).^2 );
+
+OTFpower = OTFo.*conj(OTFo);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = FiSMao.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Object Power determination
+Aobj = OBJpara(1);
+Bobj = OBJpara(2);
+Ro(wo+1,wo+1) = 1;
+OBJpower = Aobj*(Ro.^Bobj);
+OBJpower = OBJpower.^2;
+
+%% Wiener Filtering
+FiSMaof = FiSMao.*(SFo.*conj(OTFo)./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+
+%% for cross-checking filtered estimate visually
+%{
+WFilter = (SFo.^2.*OTFpower./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+FiSMao1 = FiSMaof.*OTFo;
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(FiSMaof(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao(wo+pp,:))),'o--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao1(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,10*WFilter(wo+pp,:),'k*--','LineWidth',2,'MarkerSize',4)
+title('WoFilter')
+legend('FiSMaof','FiSMao','FiSMao1')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/TIRFexpt/WoFilterSideLobe.m b/TIRFexpt/WoFilterSideLobe.m
new file mode 100644
index 0000000..fb717c5
--- /dev/null
+++ b/TIRFexpt/WoFilterSideLobe.m
@@ -0,0 +1,54 @@
+function [FiSMaof,NoisePower] = WoFilterSideLobe(FiSMao,OTFo,co,OBJsideP,SFo)
+% Aim: Wiener Filtering the off-center frequency component
+% INPUT VARIABLES
+% FiSMao: noisy off-center frequency component
+% OTFo: system OTF
+% co: Wiener filter constant [=1, for minimum RMS estimate]
+% OBJsideP: avg. object spectrum
+% SFo: modulation factor
+% OUTPUT VARIABLES
+% FiSMaof: Wiener Filtered estimate of FiSMao
+% NoisePower: avg. noise power in FiSMao
+
+w = size(FiSMao,1);
+wo = w/2;
+x = linspace(0,w-1,w);
+y = linspace(0,w-1,w);
+[X,Y] = meshgrid(x,y);
+Ro = sqrt( (X-wo).^2 + (Y-wo).^2 );
+
+OTFpower = OTFo.*conj(OTFo);
+
+% OTF cut-off frequency
+Kotf = OTFedgeF(OTFo);
+
+% frequency beyond which NoisePower estimate to be computed
+NoiseFreq = Kotf + 20;
+
+% NoisePower determination
+Zo = Ro>NoiseFreq;
+nNoise = FiSMao.*Zo;
+NoisePower = sum(sum( nNoise.*conj(nNoise) ))./sum(sum(Zo));
+
+% Object Power determination
+OBJpower = OBJsideP.^2;
+
+%% Wiener Filtering
+FiSMaof = FiSMao.*(SFo.*conj(OTFo)./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+
+%% for cross-checking filtered estimate visually
+%{
+WFilter = (SFo.^2.*OTFpower./NoisePower)./((SFo.^2).*OTFpower./NoisePower + co./OBJpower);
+FiSMao1 = FiSMaof.*OTFo;
+pp = 3;
+figure;
+hold on
+plot(x-wo,log(abs(FiSMaof(wo+pp,:))),'go--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao(wo+pp,:))),'o--','LineWidth',2,'MarkerSize',6)
+plot(x-wo,log(abs(FiSMao1(wo+pp,:))),'r*--','LineWidth',2,'MarkerSize',4)
+plot(x-wo,10*WFilter(wo+pp,:),'k*--','LineWidth',2,'MarkerSize',4)
+title('WoFilter')
+legend('FiSMaof','FiSMao','FiSMao1')
+grid on
+box on
+%}
\ No newline at end of file
diff --git a/TIRFexpt/sim01z4.tif b/TIRFexpt/sim01z4.tif
new file mode 100644
index 0000000..e8fa42d
Binary files /dev/null and b/TIRFexpt/sim01z4.tif differ