

%%
figure('Position',[500 50 740 500])
subplot(3,6,3)
fratio=[];
lfreqs=[1 2 4 8 16 32];

lfreqs=[1 2 4 7 10 15 20];


t=0:.001:1;

%IRF

A=2
B=.05
Yg = gampdf(t,A,B);
Yg=Yg/sum(Yg);

plot(t,Yg,'Color',[0 .5 1],'Linewidth',2)
xlabel('Time (s)')
set(gca,'YTick',[])
box off
xlim([0 0.5])
set(gca,'XTick',[0:.1:.5],'TickDir','out')
set(gca,'FontSize',9)
% set(gca,'PlotBoxAspectRatio',[2 3.5 2])
title('impulse response function','FontSize',9)
axis off

%signals

for i=1:length(lfreqs)

    tf=lfreqs(i);


    %squarewaves
    y=sin(2*pi*tf*t-(pi/2));
    y=(y+1)/2;
    y=round(y);

    subplot(3*length(lfreqs),3,3*(i-1)+1)
    plot(t,y,'k')
    set(gca,'YTick',[])
    set(gca,'FontSize',9)
    set(gca,'XTick',[0 0.5 1])
    if i==1
        title('input signal {\it S_f}','FontSize',9)
    end
    axis off

    %filtering
    c=conv(y,Yg);
    c=c/max(c(:));
    fratio=[fratio 1000/sum(c)];

    subplot(3*length(lfreqs),2,2*(i-1)+2)
    plot([t t(2:end)+t(end)],c,'Color',[0 .5 1])
    xlim([0 1.5])
    set(gca,'FontSize',9)
    set(gca,'YTick',[])
    set(gca,'XTick',[0 0.5 1 1.5])
    axis off
    if i==1
        title('response signal {\it E_f}','FontSize',9)
    end
end

% display timescale
i=length(lfreqs)+1;
subplot(3*length(lfreqs),3,3*(i-1)+1)
xlim([0 1])
set(gca,'FontSize',9)
set(gca,'XTick',[0 0.5 1 ],'TickDir','out')
set(gca,'YTick',[])
xlabel('Time (s)')

subplot(3*length(lfreqs),6,6*(i-1)+3)
xlim([0 .5])
set(gca,'FontSize',9)
set(gca,'XTick',[0:.1:.5 ],'TickDir','out','XTickLabel',{'0',[],[],[],[],'0.5'},'XTickLabelRotation',0 )
set(gca,'YTick',[])
xlabel('Time (s)')

subplot(3*length(lfreqs),2,2*(i-1)+2)
xlim([0 1.5])
set(gca,'FontSize',9)
set(gca,'XTick',[0 0.5 1 1.5],'TickDir','out')
set(gca,'YTick',[])
xlabel('Time (s)')

% balance point
subplot(2,2,3)
plot(lfreqs,-log(fratio),'o-','MarkerSize',12,'Color',[0 .5 1],'MarkerFaceColor',[0 .5 1])
set(gca,'XScale','log')
xlim([.7 30])
ylim([-0.8 0.1])
set(gca,'YGrid','on')
set(gca,'XTick',lfreqs)
xlabel('Temporal frequency (Hz)')
ylabel('logBP')
set(gca,'FontSize',9)

title('balance point','FontSize',9)


