% companion matlab script to lowpassexamp.m to illustrate an
% alternative to fftshift.

clear all
close all
fc=0.1;
wc=2*pi*fc;
rate=100;
fmax=rate/2;
T=1;
N=rate*T;
dt=1/rate;
df=2*fmax/N;
f=0:df:fmax;
t=0:dt:T-dt;

% now we're only doing the positive side of the response
halfphif=wc./(i*2*pi*f+wc);

subplot(2,1,1)
plot(f,abs(halfphif))
title(['positive half 1-pole ',num2str(fc),' Hz lowpass filter'])
ylabel('amplitude response')
subplot(2,1,2)
plot(f,180*angle(halfphif)/pi)
ylabel('phase response (degrees)')
xlabel('frequency (Hz)')

disp('positive half 1-pole lowpass filter');
disp(' ')
disp('press enter to continue');
disp(' ')
pause;


%orient tall
%print -dps lowpassexamp2.halfphif.ps

% flip the positive half, take the complex conjugate, and paste it
% on the right of the response vector.  The result should be Hermitian.
% But be carefule to put the DC point at the nyquist not fmax
phif=[halfphif,conj(fliplr(halfphif(:,2:N/2)))];

% plot the 2-sided xfer function 
figure
subplot(2,1,1)
plot([f,fmax+f(:,1:N/2-1)],abs(phif))
title(['two-sided 1-pole ',num2str(fc),' Hz lowpass filter'])
ylabel('amplitude response')
subplot(2,1,2)
plot([f,fmax+f(:,1:N/2-1)],180*angle(phif)/pi)
ylabel('phase response (degrees)')
xlabel('frequency (Hz)')

disp('two sided 1-pole lowpass filter');
disp(' ')
disp('press enter to continue');
disp(' ')
pause;

%orient tall
%print -dps lowpassexamp2.phif.ps

phit=ifft(phif);
figure
subplot(2,1,1)
plot(t,real(phit))
title('example two time domain transfer function of same filter')
ylabel('real part')
subplot(2,1,2)
plot(t,imag(phit))
ylabel('imaginary part')
xlabel('time (seconds)')

disp('time domain response of 1-pole lowpass filter');
disp(' ')
disp('press enter to continue');
disp(' ')

%orient tall
%print -dps lowpassexamp2.phit.ps
