Hi,
My goal is to normalize out a linear phase shift from the Fourier transform of a
rectangle with settable delay and width.
I have a strange phenomena that I have not been able to track down. I want to
take the fft of a rectangular pulse (with settable delay (tah) and width (pw))
and then look at the magnitude and phase response versus time. Mathematically
taking my representation of the rect (i.e. not strictly centered at DC) shows
that there is a linear phase shift.
exp(-jw*(pw/2 + tah)) * pw*sinc(w*(pw/2))
Mathematically there is a linear phase shift, as shown above) represented by:
exp(-jw*(pw/2 + tah))
which I normalize out with:
exp(+jw*(pw/2 + tah))
The problem is that I always have a residual phase shift that I cannot remove.
To completely flatten out the linear phase shift I need to apply:
exp(+jw*(pw/2 + tah - sr/2))
where sr is the sample rate. This works all the time no matter what the pulse
width or delay is.
I have not been able to track down this phenomena of the unknown delay of sr/2.
Below is some sample code that illustrates this. Any ideas anyone has would be
greatly appreciated.
sr = 1e-12; % Sample rate
N = 10; % Number of points in pulse width
on_time = ones(1,N); % Rect integer width
DC = .01; % Duty cycle
O = ceil(abs((N/DC) - N)); % Number of zeros in the off state
off_time = zeros(1,O); % Array of zeros in the off state
P = 0; % Number of delay zeros at beginning of pulse
delay_time = zeros(1,P); % Array of delay zeros at the beginning of pulse
amplitude = 1; % Amplitude of rect
over_sample = 1e1; % Oversample for single rect fft
unknown_delay = 0; % ****unknown delay:sr/2*****
% Scaling (Make sure to change labels when modifying these values)
t_scale = 1e-12;
f_scale = 1e9;
null_scale = 1e9;
prf_scale = 1e6;
x = amplitude.*[on_time off_time]; % Rect in time domain
if P > 0 % Initial delay?
x = [delay_time x];
end
x_size = length(x);
tah = P*sr; % Delay time
pw = N*sr; % Pulse width
null_points = 1/pw; % Sinc null points
t = 0:1:x_size-1; % Create time points
t = t*sr; % Scale time axis
% Plot time domain response
figure('color','white');
plot(t/t_scale,x);
grid on;
xlabel('Time (ps)');
ylabel('Amplitude (V)');
title(sprintf('Single DC pulse (PW = %g ps)',pw/t_scale));
% Calculate FFT (Strictly for computing phase response)
fft_points = 2^(nextpow2(over_sample*x_size)); % Power of 2. Creates symmetric
spectrum
X = fft(x,fft_points);
[PX,px] = parseval(x); % Check Parseval's theorem
unique_points = ceil((fft_points+1)/2);
X = X(1:unique_points); % Take half spectrum
X = X/x_size;
f = (0:unique_points-1)/(sr*fft_points);
% Plot frequency domain magnitude spectrum
figure('color','white');
subplot(2,1,1);
plot(f/f_scale,real(X));
grid on;
xlabel('Frequency (GHz)');
ylabel('Real');
title(sprintf('"Real" response in frequency domain(Null points are %g
GHz)',null_points/null_scale));
subplot(2,1,2);
plot(f/f_scale,imag(X));
grid on;
xlabel('Frequency (GHz)');
ylabel('Imaginary');
title(sprintf('"Imaginary" response in frequency domain(Null points are %g
GHz)',null_points/null_scale));
% Plot frequency domain phase response
figure('color','white');
plot(f/f_scale,angle(X));
grid on;
xlabel('Frequency (GHz)');
ylabel('Phase (radians)');
title(sprintf('Phase response in frequency domain(Null points are %g
GHz)',null_points/null_scale));
% Unwrap linear phase
phase_lin_factor = exp(+1j*(2*pi*f)*(((pw/2)+tah - unknown_delay))); % Angular
frequency units are radians/s
X1 = phase_lin_factor.*X;
% Plot frequency domain magnitude spectrum with linearization
figure('color','white');
subplot(2,1,1);
plot(f/f_scale,real(X1));
grid on;
xlabel('Frequency (GHz)');
ylabel('Real');
title(sprintf('"Real" response in frequency domain(Null points are %g
GHz)',null_points/null_scale));
subplot(2,1,2);
plot(f/f_scale,imag(X1));
grid on;
xlabel('Frequency (GHz)');
ylabel('Imaginary');
title(sprintf('"Imaginary" response in frequency domain(Null points are %g
GHz)',null_points/null_scale));
% Plot frequency domain phase response with linearization
figure('color','white');
plot(f/f_scale,angle(X1));
grid on;
xlabel('Frequency (GHz)');
ylabel('Phase (radians)');
title(sprintf('Phase response (Unwrap) in frequency domain(Null points are
%g GHz)',null_points/null_scale'));
