(* <pre>
Masking Model for 98ecvp  *)


(* General *)
Ds[x_]:= Dimensions[x]

Sqr[x_] := x*x

Energy[seq_] := Apply[Plus, Sqr[Flatten[seq]]]

(* generate mask image in the space domain *)

genMaskIm[f_] := 1.0+amp0 N[Sin[f w0Tab ] ]

genMaskIm[f1_,f2_] := 1.0+amp0 N[ Sin[f1 w0Tab ] + Sin[f2 w0Tab ] ]
  
(*
amp0 = 0.15 ; freq0 = 1. ; cols = 12 ;
w0 = N[2. Pi freq0/cols] ;
w0Tab =  w0 Table[ ix,{ix,0,cols-1}] ;
(* 12 points allow a frequency of 3 cycles per image in sin phase *)
x1 = genMaskIm[1.] ;
{1., 1.075, 1.1299, 1.15, 1.1299, 1.075, 1., 0.925, 0.870096,
 0.85, 0.870096, 0.925}

x3 = genMaskIm[3.] ;
{1., 1.15, 1., 0.85, 1., 1.15, 1., 0.85, 1., 1.15, 1., 0.85 } ;
(x3-1.).(x3-1.) + (x1-1.).(x1-1.) ;  0.27;

   cols = 12; freq0 = 1.; amp0 = 0.15 ;
 x13 = genMaskIm[1.,3.] ;
 {1., 1.225, 1.1299,   1., 1.1299,   1.225,
  1., 0.775, 0.870096, 1., 0.870096, 0.775}
 (x13-1.).(x13-1.) = 0.269998 ;
 ListPlot[x13, Prolog ->{PointSize[0.02]}]
*)
   
(* generate signal image in the space domain *)

genSigIm[f_] := sigAmp N[Sin[f  w0Tab ] ]

(* generate sequence of 1 dimensional grating images in the
   spatial frequency domain.
   fTab0 = spatial frequency image for beginning and end
   fTab1 = spatial frequency image for mask region 
   lengths = {frames in start, frames in mask, frames at end}
   cols is the number of pixels in each image.
   w0 is the spatial frequency in radians per pixel.
   amp0 is the peak contrast of the grating.
   *)

genMaskSeq := (
        fSeq = Join[
         Table[ fTab0, {lengths[[1]]}],
         Table[ fTab1, {lengths[[2]]}],
         Table[ fTab0, {lengths[[3]]}] ]);

(* Filtering *)

(* temporal *)

(* convert exponential decay filter cutoff f0 in frequency domain
   Amp[f] = f0/Sqrt[f f + f0 f0]
   to digital attenuation factor *)

expCut2DigAtten[f0_, fs_] := Exp[N[-2.*Pi*f0/fs] ] ;

(* low pass filter *)
tLP[seq_,a_,y0_] :=Module[{i,y1,out},
 y1=y0;
 out = seq;
 For[i=1, i<1+Length[seq], i++,
    y1=(1.-a)*seq[[i]]+a*y1;
    out[[i]] =  y1
 ];
 out ];

(* high pass filter; successive differences scaled by a *)
tHP[seq_,a_,x0_] :=Module[{i,x1,x2,out},
 x1=x0;
 out = seq;
 For[i=1, i<1+Length[seq], i++,
    out[[i]] = a*((x2=seq[[i]])-x1);
    x1=x2;
 ];
 out ];

(* test of tLP
seq = {1.,0.,0.}; a = 0.3; y0 = 0.;
tLP[seq,a,y0] ; {0.7, 0.21, 0.063} ;
*)

seta := (
  ac = expCut2DigAtten[ftc, fs] ;
  ah = expCut2DigAtten[fth, fs] ;
  ap = expCut2DigAtten[ftp, fs] ; )

(* spatial *)

(* nx must be even *)
sFilterG1[nx_, sx_] := 
  Chop[Exp[-(Join[Table[x,{x,0,nx/2  }],
                  Table[x,{x,1,nx/2-1}]-nx/2]/sx)^2] ] ;

setFilts := (
  filtc = sFilterG1[cols, fsc* cols/fres ];
  filth = sFilterG1[cols, fsh* cols/fres ];
  filtm = sFilterG1[cols, fsm* cols/fres ];
  filtq = sFilterG1[cols, fsq* cols/fres ];
  )

(* spatiotemporal *)

(* frequency sequence to low pass spatial *)
f2LPs[fSeq_, filt_, a_, y0_] :=
  tLP[Map[Re[Fourier[filt*#] ]&,fSeq], a, y0]

(* generate model responses *)

(* models use global variables to compute parvo and magno outputs:
   fSeq = input sequence in frequency domain
   fSeq0 = input sequence initializing image in frequency domain
   filtc, filth, filtm = spatial filters
   ac, ah, ap = temporal filter parameters
   cone0, coneblur, cone, horiz, bipolar, parvo= intermediate outputs
 *)

(* no masking, no temporal filtering in magno *)
model :=(
  cone0 = Re[Fourier[filtc*fSeq0] ] ;
  coneblur = Map[(Re[Fourier[filtc*#] ])&,fSeq] ;
  cone = tLP[coneblur, ac, cone0 ] ;
 horiz = f2LPs[fSeq, filth, ah, cone0 ]  ;
 bipolar = N[ cone / horiz - 1. ] ;
 parvo = N[tLP[bipolar, ap, zeros] ] ;
 magno = Map[(Re[Fourier[filtm*InverseFourier[#]]])&,bipolar] ;
 )

(* no masking, temporal filtering in magno is LP,HP *)
model2 :=(
  zeros = Table[0.0,{cols}] ;
  cone0 = Re[Fourier[filtc*fSeq0] ] ;
  coneblur = Map[(Re[Fourier[filtc*#] ])&,fSeq] ;
  cone = tLP[coneblur, ac, cone0 ] ;
 horiz = f2LPs[fSeq, filth, ah, cone0 ]  ;
 bipolar = N[ cone / horiz - 1. ] ;
 parvo = N[tLP[bipolar, ap, zeros] ] ;
 (*
 magno = Map[(Re[Fourier[filtm*InverseFourier[#]]])&,bipolar] ;
 magno = N[tLP[magno, ap, zeros] ] ;
 *)
 magno = Map[(Re[Fourier[filtm*InverseFourier[#]]])&,parvo] ;
 magno = N[tHP[magno, 1., zeros] ] ;
 )

(* no masking, temporal filtering in magno is LP,LP,HP,HP *)
(*
model3 :=(
  cone0 = Re[Fourier[filtc*fSeq0] ] ;
  coneblur = Map[(Re[Fourier[filtc*#] ])&,fSeq] ;
  cone = tLP[coneblur, ac, cone0 ] ;
 horiz = f2LPs[fSeq, filth, ah, cone0 ]  ;
 bipolar = N[ cone / horiz - 1. ] ;
 parvo = N[tLP[bipolar, ap, zeros] ] ;
 magno = Map[(Re[Fourier[filtm*InverseFourier[#]]])&,bipolar] ;
 magno = N[tLP[magno, ap, zeros] ] ;
 magno = N[tLP[magno, ap, zeros] ] ;
 magno = N[tHP[magno, 1., zeros] ] ;
 magno = N[tHP[magno, 1., zeros] ] ;
 )
*)
maskF[seq_,filt_,sig_] := seq/Sqrt[1.+
          Map[(Re[Fourier[filt*InverseFourier[#]]])&,Sqr[seq/sig] ] ];

mask0[seq_,sig_] := seq/Sqrt[1.+ Sqr[seq/sig] ] ;

maskT[seq_,sig_,a_] := N[seq/Sqrt[1.+ tLP[Sqr[seq/sig],a,zeros] ] ] ;

maskTF[seq_,filt_,sig_,a_] := N[seq/Sqrt[1.+ tLP[
          Map[(Re[Fourier[filt*InverseFourier[#]]])&,Sqr[seq/sig] ],
          a,zeros] ] ];


maskModel2 :=(
 model2;
 magnom = maskF[magno,filtq,sigm] ;
 parvom = maskF[parvo,filtq,sigp] ;
 )

(*
maskModel3 :=(
 model3;
 magnom = maskF[magno,filtq,sigm] ;
 parvom = maskF[parvo,filtq,sigp] ;
 )
*)

maskModel4 :=(
 model2;
 magnom = maskTF[magno,filtq,sigm,ap] ;
 parvom = maskT[parvo,sigp,ap] ;
 )

maskModel4 :=(
 model2;
 magnom = maskTF[magno,filtq,sigm,ap] ;
 parvom = mask0[parvo,sigp] ;
 )

(*-----------------------------------------------*)


(* </pre> *)