P=PLExp1_params ; [Fx,Fy]=fftmesh(P.size(1)); Fx=fftshift(Fx)/P.deg_pix;Fy=fftshift(Fy)/P.deg_pix;
o=[-45:2.5:+45];orient=o*(pi/180);sp_freq=[0:18]'/3;method='*cubic';idx=[10 8 7];

%%% It's better to express the spatial frequencies on a logarithmic scale of "octaves":
%log_sp_freq=[-1:0.25:+3]';sp_freq=2.^log_sp_freq;idx=[11 10 9];

%%% Generate a stimulus, calculate its spectrum, convert to polar coordinates, and plot
x=make_PLE_stim(P,1,1,2)-P.colors.gray;pwr=fftshift(abs(fft2(x/128)));rep=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);subplot(2,1,1);imagesc(orient*180/pi,sp_freq,rep,[0 80]);colorbar;subplot(2,1,2);plot(o,rep(idx,:),'.-');title([method ';  sp. freq = ' num2str(sp_freq(idx(3)))]);axis([min(o) max(o) 0 80]);grid on;

%%% Several individual instances
for k=1:30 foo=2-(mod(k-1,6)<3);x=make_PLE_stim(P,foo,1,1)-128;pwr=fftshift(abs(fft2(x/128)));rep=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);subplot(5,6,k);plot(o,rep(7,:),'b-');axis([min(o) max(o) 0 80]);set(gca,'xtick',[-10 +10],'ytick',[0:20:80]);grid on;end

%%% Collapse the spectrum to a two-element vector and produce scatterplots
neg_idx=find(orient<0);pos_idx=find(orient>0);x=make_PLE_stim(P,1,1,2)-P.colors.gray;pwr=fftshift(abs(fft2(x/128)));rep=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);L=mean(rep(7,neg_idx));R=mean(rep(7,pos_idx));[L R]
LR=zeros(100,2);for k=1:100 x=make_PLE_stim(P,1,1,1)-P.colors.gray;pwr=fftshift(abs(fft2(x/128)));rep=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);LR(k,1)=mean(rep(7,neg_idx));LR(k,2)=mean(rep(7,pos_idx));end
LR_1_1_1=LR;           % and similarly for LR_1_1_2, etc.
plot(LR(:,1),LR(:,2),'.');axis([0 60 0 60]);refline(1,0);xlabel('Avg power at negative orient');ylabel('AvgPwr at pos orient');

%%% Discrimination
plot(LR_1_3_1(:,1),LR_1_3_1(:,2),'b.',LR_2_3_1(:,1),LR_2_3_1(:,2),'r.');axis([0 60 0 60]);refline(1,0);xlabel('Avg power at negative orient');ylabel('AvgPwr at pos orient');

plot(LR_1_1_1(:,1),LR_1_1_1(:,2),'b.',LR_2_1_1(:,1),LR_2_1_1(:,2),'r.');axis([0 60 0 60]);refline(1,0);xlabel('Avg power at negative orient');ylabel('AvgPwr at pos orient');
hold on;plot(LR_1_2_1(:,1),LR_1_2_1(:,2),'b.',LR_2_2_1(:,1),LR_2_2_1(:,2),'r.');hold off
hold on;plot(LR_1_3_1(:,1),LR_1_3_1(:,2),'b.',LR_2_3_1(:,1),LR_2_3_1(:,2),'r.');hold off


%%% Illustrate the discrimination task, noisy representations
x=make_PLE_stim(P,1,1,1)-128;pwr=fftshift(abs(fft2(x/128)));rep1=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);x=make_PLE_stim(P,2,1,1)-128;pwr=fftshift(abs(fft2(x/128)));rep2=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);
plot(o,rep1(7,:),'r.-',o,rep2(7,:),'b.:');axis([min(o) max(o) 0 80]);grid on;xlabel('Orientation [degrees]');ylabel('Spectral power');title('Noise -15 deg ; Sp.freq 2 cyc/deg');legend('Gabor -10, C_p=0.23','Gabor +10, C_p=0.23');

for k=1:12 x=make_PLE_stim(P,1,1,1)-128;pwr=fftshift(abs(fft2(x/128)));rep1=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);x=make_PLE_stim(P,2,1,1)-128;pwr=fftshift(abs(fft2(x/128)));rep2=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);subplot(3,4,k);plot(o,rep1(7,:),'r.-',o,rep2(7,:),'b.:');axis([min(o) max(o) 0 80]);end

%%% Average over 100 Monte Carlo trials
foo=zeros(19,37,100);for k=1:100 x=make_PLE_stim(P,1,1,2)-P.colors.gray;pwr=fftshift(abs(fft2(x/128)));rep=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);foo(:,:,k)=rep;end
rep=squeeze(mean(foo,3));rep_1_1_2=rep;
subplot(2,1,1);imagesc(orient*180/pi,sp_freq,rep,[0 80]);colorbar;subplot(2,1,2);plot(o,rep(idx,:),'.-');title([method ';  sp. freq = ' num2str(sp_freq(idx(3)))]);axis([min(o) max(o) 0 80]);grid on;
%.. analogously for rep_1_1_1, rep_2_1_1, etc.

%%%  Check for symmetry
rep=rep_1_1_1-fliplr(rep_2_1_2);subplot(2,1,1);imagesc(orient*180/pi,sp_freq,rep,[-5 5]);colorbar;subplot(2,1,2);plot(o,rep(idx,:),'.-');title([method ';  sp. freq = ' num2str(sp_freq(idx(3)))]);axis([min(o) max(o) -5 5]);grid on;

%%%  Plot a congruent and an incongruent stimulus (maximal Gabor contrast)
plot(o,rep_1_1_1(7,:),'r.-',o,rep_1_1_2(7,:),'b.:');axis([min(o) max(o) 0 80]);grid on;xlabel('Orientation [degrees]');ylabel('Spectral power');title('Gabor contrast 0.23; Sp.freq 2 cyc/deg');legend('congruent','incongruent');
plot(o,rep_1_2_1(7,:),'r.-',o,rep_1_2_2(7,:),'b.:');axis([min(o) max(o) 0 80]);grid on;xlabel('Orientation [degrees]');ylabel('Spectral power');title('Gabor contrast 0.15; Sp.freq 2 cyc/deg');legend('congruent','incongruent');
plot(o,rep_1_3_1(7,:),'r.-',o,rep_1_3_2(7,:),'b.:');axis([min(o) max(o) 0 80]);grid on;xlabel('Orientation [degrees]');ylabel('Spectral power');title('Gabor contrast 0.10; Sp.freq 2 cyc/deg');legend('congruent','incongruent');

plot(o,rep_1_1_1(7,:),'r.-',o,rep_2_1_1(7,:),'b.:');axis([min(o) max(o) 0 80]);grid on;xlabel('Orientation [degrees]');ylabel('Spectral power');title('Noise -15 deg ; Sp.freq 2 cyc/deg');legend('Gabor -10, C_p=0.23','Gabor +10, C_p=0.23');
plot(o,rep_1_2_1(7,:),'r.-',o,rep_2_2_1(7,:),'b.:');axis([min(o) max(o) 0 80]);grid on;xlabel('Orientation [degrees]');ylabel('Spectral power');title('Noise -15 deg ; Sp.freq 2 cyc/deg');legend('Gabor -10, C_p=0.15','Gabor +10, C_p=0.15');
plot(o,rep_1_3_1(7,:),'r.-',o,rep_2_3_1(7,:),'b.:');axis([min(o) max(o) 0 80]);grid on;xlabel('Orientation [degrees]');ylabel('Spectral power');title('Noise -15 deg ; Sp.freq 2 cyc/deg');legend('Gabor -10, C_p=0.10','Gabor +10, C_p=0.10');

%%%  Subtract two contrast levels
rep=rep_1_1_2-rep_1_3_2;subplot(2,1,1);imagesc(orient*180/pi,sp_freq,rep,[0 80]);colorbar;subplot(2,1,2);plot(o,rep(idx,:),'.-');title([method ';  sp. freq = ' num2str(sp_freq(idx(3)))]);axis([min(o) max(o) 0 80]);grid on;


%%%  Estimate the spectrum of the NOISE alone
Cp0=length(P.Gabor_Cpeak)+1;P.Gabor_Cpeak(Cp0)=0;P.Gabor(:,:,:,Cp0)=zeros([P.size,length(P.Gabor_orient)]);
foo=zeros(19,37,100);for k=1:100 x=make_PLE_stim(P,1,Cp0,1)-P.colors.gray;pwr=fftshift(abs(fft2(x/128)));rep=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);foo(:,:,k)=rep;end
rep=squeeze(mean(foo,3));rep_noise_1=rep;
subplot(2,1,1);imagesc(orient*180/pi,sp_freq,rep,[0 80]);colorbar;subplot(2,1,2);plot(o,rep(idx,:),'.-');title([method ';  sp. freq = ' num2str(sp_freq(idx(3)))]);axis([min(o) max(o) 0 80]);grid on;

%%%  Check the noise for symmetry
rep=rep_noise_1-fliplr(rep_noise_2);subplot(2,1,1);imagesc(orient*180/pi,sp_freq,rep,[-5 5]);colorbar;subplot(2,1,2);plot(o,rep(idx,:),'.-');title([method ';  sp. freq = ' num2str(sp_freq(idx(3)))]);axis([min(o) max(o) -5 5]);grid on;
plot(o,rep_noise_1(7,:),'r.-',o,rep_noise_2(7,:),'b.:');axis([min(o) max(o) 0 80]);grid on;xlabel('Orientation [degrees]');ylabel('Spectral power');title('Noise alone;  Spatial frequency 2 cycles/degree');legend('noise 1','noise 2');

%%%  Individual noise spectra
for k=1:20 x=make_PLE_stim(P,1,Cp0,1)-128;pwr=fftshift(abs(fft2(x/128)));rep=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);subplot(4,5,k);plot(o,rep(7,:),'b.-');axis([min(o) max(o) 0 60]);set(gca,'xtick',[-15 0 +15],'ytick',[0:20:60]);grid on;end

%%%  Partial out the noise from the composite stimuli
rep=rep_1_1_1-rep_noise_1;subplot(2,1,1);imagesc(orient*180/pi,sp_freq,rep,[0 80]);colorbar;subplot(2,1,2);plot(o,rep(idx,:),'.-');title([method ';  sp. freq = ' num2str(sp_freq(idx(3)))]);axis([min(o) max(o) -20 60]);grid on;


%%%  Gabor alone
noise_Cpeak=P.noise_Cpeak;P.noise_Cpeak=0;x=make_PLE_stim(P,1,1,1)-P.colors.gray;P.noise_Cpeak=noise_Cpeak;pwr=fftshift(abs(fft2(x/128)));rep=cart2pol2(Fx,Fy,pwr,orient,sp_freq,method);subplot(2,1,1);imagesc(orient*180/pi,sp_freq,rep,[0 80]);colorbar;subplot(2,1,2);plot(o,rep(idx,:),'.-');title([method ';  sp. freq = ' num2str(sp_freq(idx(3)))]);axis([min(o) max(o) 0 80]);grid on;
rep_Gabor_1_1=rep;   % ... and so on for Gabor_1_2, etc.
figure;plot(o,[rep_Gabor_1_1(7,:);rep_Gabor_1_2(7,:);rep_Gabor_1_3(7,:)]','.-');axis([min(o) max(o) 0 80]);grid on;xlabel('Orientation [degrees]');ylabel('Spectral power');title('Gabors alone;  Spatial frequency 2 cycles/degree');legend(['C_p=',num2str(P.Gabor_Cpeak(1))],['C_p=',num2str(P.Gabor_Cpeak(2))],['C_p=',num2str(P.Gabor_Cpeak(3))]);
foo=[rep_Gabor_1_1(7,:);rep_Gabor_1_2(7,:);rep_Gabor_1_3(7,:)]';
[foo(15,:);P.Gabor_Cpeak(1:3);foo(15,:)./P.Gabor_Cpeak(1:3)]    % peak values, orient=-10 deg



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%   OBSOLETE   %%%%%%%%%%%%%%%%%%%%%%


clear all ; PLExp1 ;   % A RETURN clause must be included at the appropriate point

% Print one sample of each stimulus configuration
k=1;for a=1:2 for b=1:2 subplot(2,2,k);img=make_PLE_stim(P,a,1,b);imagesc(img);axis('image');colormap gray;set(gca,'xtick',[],'ytick',[]);k=k+1;end;end

% Estimate the spectrum via 100 Monte Carlo samples in each condition
pStim=zeros(64,64,2,2);for a=1:2 for b=1:2 p=zeros(64,64,100);for k=1:100 img=(make_PLE_stim(P,a,1,b)-128)./256;p(:,:,k)=log(max(4,fftshift(abs(fft2(img)))));end;pStim(:,:,a,b)=squeeze(mean(p,3));end;end

% plot in color using IMAGESC
k=1;for a=1:2 for b=1:2 subplot(2,2,k);imagesc(pStim(:,:,a,b));colorbar;axis('image');title(sprintf('Gabor %d  ; noise %d',P.Gabor_orient(a)*180/pi,round(P.noise_orient(b)*180/pi)));k=k+1;colormap jet;set(gca,'xtick',[],'ytick',[]);end;end

% plot in white-and-black for printing
k=1;for a=1:2 for b=1:2 subplot(2,2,k);imagesc(-pStim(:,:,a,b));colorbar;axis('image');title(sprintf('Gabor %d  ; noise %d',P.Gabor_orient(a)*180/pi,round(P.noise_orient(b)*180/pi)));k=k+1;colormap gray;set(gca,'xtick',[],'ytick',[]);end;end

% contour plot in color   (Note that the data are flipped upside down because CONTOURF puts Y=0 at bottom)
k=1;for a=1:2 for b=1:2 subplot(2,2,k);contourf(pStim(64:-1:1,:,a,b));colorbar;axis('image');title(sprintf('Gabor %d  ; noise %d',P.Gabor_orient(a)*180/pi,round(P.noise_orient(b)*180/pi)));k=k+1;colormap jet;set(gca,'xtick',[],'ytick',[]);end;end

% contour plot in BW
k=1;for a=1:2 for b=1:2 subplot(2,2,k);contourf(-pStim(64:-1:1,:,a,b));colorbar;axis('image');title(sprintf('Gabor %d  ; noise %d',P.Gabor_orient(a)*180/pi,round(P.noise_orient(b)*180/pi)));k=k+1;colormap gray;set(gca,'xtick',[],'ytick',[]);end;end
