5 files changed, 162 insertions, 0 deletions
diff --git a/Measurements/A_f.m b/Measurements/A_f.m
new file mode 100644
index 0000000..3a793aa
--- /dev/null
+++ b/Measurements/A_f.m
@@ -0,0 +1,30 @@
+% A_f.m
+%
+% Takes "scrambled Fourier" measurements.
+%
+% Usage: b = A_f(x, OMEGA, P)
+%
+% x - N vector
+%
+% b - K vector = [real part; imag part]
+%
+% OMEGA - K/2 vector denoting which Fourier coefficients to use
+%         (the real and imag parts of each freq are kept).
+%
+% P - Permutation to apply to the input vector.  Fourier coeffs of
+%     x(P) are calculated.
+%     Default = 1:N (no scrambling).
+%
+% Written by: Justin Romberg, Caltech
+% Created: October 2005
+% Email: jrom@acm.caltech.edu
+%
+
+function b = A_f(x, OMEGA, P)
+
+N = length(x);
+if (nargin < 3), P = 1:N;  end
+
+fx = 1/sqrt(N)*fft(x(P));
+b = [sqrt(2)*real(fx(OMEGA)); sqrt(2)*imag(fx(OMEGA))];
+
diff --git a/Measurements/A_fhp.m b/Measurements/A_fhp.m
new file mode 100644
index 0000000..a7aef36
--- /dev/null
+++ b/Measurements/A_fhp.m
@@ -0,0 +1,25 @@
+% A_fhp.m
+%
+% Takes measurements in the upper half-plane of the 2D Fourier transform.
+%
+% Usage: b = A_fhp(x, OMEGA)
+%
+% x - N vector
+%
+% b - K vector = [mean; real part(OMEGA); imag part(OMEGA)]
+%
+% OMEGA - K/2-1 vector denoting which Fourier coefficients to use
+%         (the real and imag parts of each freq are kept).
+%
+% Written by: Justin Romberg, Caltech
+% Created: October 2005
+% Email: jrom@acm.caltech.edu
+%
+
+function y = A_fhp(x, OMEGA)
+
+n = round(sqrt(length(x)));
+
+yc = 1/n*fft2(reshape(x,n,n));
+y = [yc(1,1); sqrt(2)*real(yc(OMEGA)); sqrt(2)*imag(yc(OMEGA))];
+
diff --git a/Measurements/At_f.m b/Measurements/At_f.m
new file mode 100644
index 0000000..9addfe3
--- /dev/null
+++ b/Measurements/At_f.m
@@ -0,0 +1,32 @@
+% At_f.m
+%
+% Adjoint for "scrambled Fourier" measurements.
+%
+% Usage: x = At_f(b, N, OMEGA, P)
+%
+% b - K vector = [real part; imag part]
+%
+% N - length of output x
+%
+% OMEGA - K/2 vector denoting which Fourier coefficients to use
+%         (the real and imag parts of each freq are kept).
+%
+% P - Permutation to apply to the input vector.  Fourier coeffs of
+%     x(P) are embedded.
+%     Default = 1:N (no scrambling).
+%
+% Written by: Justin Romberg, Caltech
+% Created: October 2005
+% Email: jrom@acm.caltech.edu
+%
+
+
+function x = At_f(b, N, OMEGA, P)
+
+if (nargin < 4),  P = 1:N;  end
+
+K = length(b);
+fx = zeros(N,1);
+fx(OMEGA) = sqrt(2)*b(1:K/2) + i*sqrt(2)*b(K/2+1:K);
+x = zeros(N,1);
+x(P) = sqrt(N)*real(ifft(fx));
diff --git a/Measurements/At_fhp.m b/Measurements/At_fhp.m
new file mode 100644
index 0000000..a7879af
--- /dev/null
+++ b/Measurements/At_fhp.m
@@ -0,0 +1,28 @@
+% At_fhp.m
+%
+% Adjoint of At_fhp (2D Fourier half plane measurements).
+%
+% Usage: x = At_fhp(b, OMEGA, n)
+%
+% b - K vector = [mean; real part(OMEGA); imag part(OMEGA)]
+%
+% OMEGA - K/2-1 vector denoting which Fourier coefficients to use
+%         (the real and imag parts of each freq are kept).
+%
+% n - Image is nxn pixels
+%
+% x - N vector
+%
+% Written by: Justin Romberg, Caltech
+% Created: October 2005
+% Email: jrom@acm.caltech.edu
+%
+
+function x = At_fhp(y, OMEGA, n)
+
+K = length(y);
+
+fx = zeros(n,n);
+fx(1,1) = y(1);
+fx(OMEGA) = sqrt(2)*(y(2:(K+1)/2) + i*y((K+3)/2:K));
+x = reshape(real(n*ifft2(fx)), n*n, 1);
diff --git a/Measurements/LineMask.m b/Measurements/LineMask.m
new file mode 100644
index 0000000..d7e985f
--- /dev/null
+++ b/Measurements/LineMask.m
@@ -0,0 +1,47 @@
+% LineMask.m
+%
+% Returns the indicator of the domain in 2D fourier space for the 
+% specified line geometry.
+% Usage :  [M,Mh,mi,mhi] = LineMask(L,N)
+%
+% Written by : Justin Romberg
+% Created : 1/26/2004
+% Revised : 12/2/2004
+
+function [M,Mh,mi,mhi] = LineMask(L,N)
+
+
+thc = linspace(0, pi-pi/L, L);
+%thc = linspace(pi/(2*L), pi-pi/(2*L), L);
+
+M = zeros(N);
+
+% full mask
+for ll = 1:L
+
+	if ((thc(ll) <= pi/4) | (thc(ll) > 3*pi/4))
+		yr = round(tan(thc(ll))*(-N/2+1:N/2-1))+N/2+1;
+    	for nn = 1:N-1
+      	M(yr(nn),nn+1) = 1;
+      end
+  else 
+		xc = round(cot(thc(ll))*(-N/2+1:N/2-1))+N/2+1;
+		for nn = 1:N-1
+			M(nn+1,xc(nn)) = 1;
+		end
+	end
+
+end
+
+
+% upper half plane mask (not including origin)
+Mh = zeros(N);
+Mh = M;
+Mh(N/2+2:N,:) = 0;
+Mh(N/2+1,N/2+1:N) = 0;
+
+
+M = ifftshift(M);
+mi = find(M);
+Mh = ifftshift(Mh);
+mhi = find(Mh);