% matlab script to demonstrate correctly mapping the negative frequency
% component to the upper half of the indeces using a basic low pass filter.

% clean my workspace
clear all
close all

% define the corner frequency as 0.1 Hz
fc=0.1;
wc=2*pi*fc;

% use a sample rate of 100 samples per second
rate=100;

% set the maximum allowed frequency (half the nyquist)
fmax=rate/2;

% the time step in seconds per sample
dt=1/rate;

% total length of the time series in seconds
T=1;

% make a time vector with an even number of points
t=0:dt:T-dt;

% the number of points
N=rate*T;

% set the frequency step in Hz
df=2*fmax/N;

% create a frequency vector
f=-fmax:df:fmax-df;

% get the frequency domain transfer function for the filter
phif=wc./(i*2*pi*f+wc);

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

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

% print fig to file
%orient tall
%print -dps lowpassexamp.phif.ps

% swap the left and right halves and inverse fft
phit=ifft(fftshift(phif));
my sum = sum(abs(imag(phit)));
if sum < 0.00000001
x=real(phit)

% plot the time domain impulse response of the low pass filter
figure
subplot(2,1,1)
plot(t,real(phit))
title(['time domain impulse response of 1-pole',num2str(fc),' Hz lowpass 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(' ')
pause;

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

% now do it again but don't swap left and right halves
phit=ifft(phif);
figure
subplot(2,1,1)
plot(t,real(phit))
title('incorrect time-domain impulse of low pass filter')
ylabel('real part')
subplot(2,1,2)
plot(t,imag(phit))
ylabel('imaginary part')
xlabel('time (seconds)')

disp('Time-domain response of without using fftshift');
disp(' ')
disp('press enter to continue');
disp(' ')

%orient tall
%print -dps lowpassexamp.wrong.ps

