%%%% This file calculates the orientation bandwidth of the Gabor stimuli

 P=PLExp1_params ;
 G11=P.Gabor(:,:,1,1);fG11=abs(fftshift(fft2(G11)));imagesc(fG11);axis square;colorbar
 G21=P.Gabor(:,:,2,1);fG21=abs(fftshift(fft2(G21)));imagesc(fG21);axis square;colorbar


 edit pol2cart2
 [Fx,Fy]=fftmesh(64);Fx=fftshift(Fx);Fy=fftshift(Fy);contour(Fx,Fy,fG11);colorbar;grid on

 or=[-45:+45];orient=or.*(pi/180); spfr=[0:0.25:4]';sp_freq=spfr.*P.deg_pix;
 pfG11=cart2pol2(Fx,Fy,fG11,orient,sp_freq);contour(or,spfr,pfG11);colorbar;grid on
 pfG21=cart2pol2(Fx,Fy,fG21,orient,sp_freq);contour(or,spfr,pfG21);colorbar;grid on

 g1=pfG11(9,:);g2=pfG21(9,:);plot(or,g1,'b-',or,g2,'r:');grid on;describe(g1)

    Mean  Std.dev       Min     Q25  Median     Q75     Max
------------------------------------------------------------
  17.864   18.634      0.01    1.13    9.38   34.82   55.93

 set(gca,'xtick',[-45:5:45],'ytick',[0:7:56]);

%%  Full height = 56; Half height = 28; 
%%  On the diagram, measure Full width at half height = FWHH = 28 degrees

 Idef=PLM_input_params;FWHH=28; FWHH*(pi/180)*Idef.W2sigma
  ans = 0.2075   % sigma in radians

 ans*180/pi
  ans = 11.8905  % sigma in degrees

 foo=normpdf(or,-10,11.89);foo=foo.*(max(g1)/max(foo));hold on;plot(or,foo,'g-');hold off

%%% Conclusion:  The Fourier transform of the Gabor stimulus has a radially
%%% Gaussian shape (as it should) centered on theta=-10 and spfreq=2 c/d.
%%% The sigma of this Gaussian is approximately 11.9 degrees = 0.2075 radians.


%%%%%%%%%%%%%%  Now, the same thing can be calculated theoretically
%  Excerpt from the documentation of the function RFIELD_SIZE.M
%  ...
%  Each Gaussian kernel in the frequency domain corresponds to a Gabor
%  receptive field in the space domain. RFIELD_SIZE calculates the
%  standard deviations [in deg.vis.ang] of the Gaussian envelopes of
%  these RFields. S_PERP is in the direction perpendicular to the bars
%  and S_PRLL parallel to them. Each sigma is inversely proportional to
%  the corresponding sigma in the frequency domain (see Table 2.2 in
%  Norma Graham's "Visual Pattern Analyzers", Oxford Univ.Press, 1989).

% For our stimuli, sigma_parallel=0.4 dva, R=2 cyc/deg
 sigma_orient=1/(2*pi * 0.4*2)
    0.1989   % radians

 sigma_orient*180/pi
   11.3986   % degrees, in excellent agreement with the 11.89 estimated via FFT



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
%%  Fourier transforms of actual stimuli

 or=[-45:+45];orient=or.*(pi/180); spfr=[0:0.25:6]';sp_freq=spfr.*P.deg_pix;

% Generate 1000 instances of PLE_stim(2,2,1) and take FFT slices through 2 cyc/deg and 4 cyc/deg
 foo2=zeros(1000,91);foo4=foo2;for k=1:1000 s=(make_PLE_stim(P,2,2,1)-P.colors.gray)./P.LUT_params.span;fs=abs(fftshift(fft2(s)));pfs=cart2pol2(Fx,Fy,fs,orient,sp_freq);foo2(k,:)=pfs(9,:);foo4(k,:)=pfs(17,:);end
 mA2=mean(foo2)';sA2=std(foo2)';mA4=mean(foo4)';sA4=std(foo4)';

% Do the same for a congruent stimulus: PLE_stim(1,2,1), intermediate Gabor contrast
 foo2=zeros(1000,91);foo4=foo2;for k=1:1000 s=(make_PLE_stim(P,1,2,1)-P.colors.gray)./P.LUT_params.span;fs=abs(fftshift(fft2(s)));pfs=cart2pol2(Fx,Fy,fs,orient,sp_freq);foo2(k,:)=pfs(9,:);foo4(k,:)=pfs(17,:);end
 mB2=mean(foo2)';sB2=std(foo2)';mB4=mean(foo4)';sB4=std(foo4)';

 subplot(2,2,1);plot(or,mB2,'k-',or,[mB2-sB2 mB2+sB2],'k--');axis([-45 45 0 50]);set(gca,'xtick',[-45:15:45]); title('Congruent stimulus');ylabel('Power at 2 cyc/deg');
 subplot(2,2,2);plot(or,mA2,'k-',or,[mA2-sA2 mA2+sA2],'k--');axis([-45 45 0 50]);set(gca,'xtick',[-45:15:45]); title('Incongruent stimulus');ylabel('Power at 2 cyc/deg');
 subplot(2,2,3);plot(or,mB4,'k-',or,[mB4-sB4 mB4+sB4],'k--');axis([-45 45 0 50]);set(gca,'xtick',[-45:15:45]); xlabel('Orientation');ylabel('Power at 4 cyc/deg');
 subplot(2,2,4);plot(or,mA4,'k-',or,[mA4-sA4 mA4+sA4],'k--');axis([-45 45 0 50]);set(gca,'xtick',[-45:15:45]); xlabel('Orientation');ylabel('Power at 4 cyc/deg');
% Page setup: [1.25 2.65 6 5.7] ;   Set line width 1.5 for the mean and 1.0 for std's
