% Jon Newman
% 2013 12 20

tau = [0.1 0.25 0.5 1 2.5 5 10];
tau = [0.16 0.8 1.6 2.5];
Fs = 100;
close all
w = 0:0.001:pi;

% Alpha
tau_n = 2*pi./(tau*Fs);
a = (tau_n./(tau_n + 1))';

% Angular frequency

% Amplitude response
a = a(:,ones(1,size(w,2)));
w = w(ones(size(a,1),1),:);

% Amplitude response
h_w = a./sqrt(1 - 2*(1-a).*cos(w) + (1-a).^2);

% Phase response
p_w = atan((-(1-a).*sin(w))./(1 - (1-a).*cos(w)));

close all 
figure('Units','centimeters','Position',[3 3 10 15]);
axes('Units','centimeters','Position',[3 6 6 4]);
semilogx(Fs*w'./pi,20*log10(h_w'))
hold on
semilogx(Fs*w'./pi,-3*ones(size(w')),'k--')
set(gca,'YTick',[-20 -15 -10 -5 -3 0])
ylim([-20 1])
xlim([0 50]);
ylabel('Amplitude (dB)','Fontsize',7,'FontName','Arial')
legend(num2str(tau'),'location','SouthWest')
set(gca,'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial','box','off')


axes('Units','centimeters','Position',[3 1 6 4]);
semilogx(Fs*w'./pi,180*p_w'./(pi/2))
hold on
semilogx(Fs*w'./pi,-180*ones(size(w')),'k--')
xlim([0 50]);
ylabel('Phase (deg.)','Fontsize',7,'FontName','Arial')
xlabel('Frequency (Hz)','Fontsize',7,'FontName','Arial')
set(gca,'TickDir','out','Fontsize',7,'TickLength',[0.01 0.01],'box','off','FontName','Arial','box','off')

saveas(gcf,'EMA_Freq-Response-vs-tau_00')
export_fig('-pdf','-transparent','EMA_Freq-Response-vs-tau_00',gcf)

