0
$\begingroup$

I am trying to recreate a function in Matlab that claimed to "use FFT to create an approximation to the continuous time Fourier transform" and was purportedly accurate for arbitrary input signals. The challenge is that I haven't been able to get it to work consistently for zero centered rectangular pulses or sinc functions--or so far anything but a Gaussian pulse. The central issue seems to be the fact that the FFT assumes the very first sample exists at time t=0. I am aware of the approximation formulas given in resources like A tutorial on THz pulse analysis: accurate retrieval of pulse arrival times, spectral energy density and absolute spectral phase and elsewhere that have two components:

  1. Scaling the FFT by deltaT (the sampling period in time)
  2. Multiplying by a complex exponential with a phase of 2pifreq*t(1)

I would greatly appreciate any insight, as all the resources I've found that either use the two step approach above or talk about ifftshift-ing the time-domain input before applying the FFT just don't seem to produce the right result. I believe I don't have any toolbox-specific functions in this script, but I can update it if so.

For the rectangular pulse, I'm expecting the phase spectrum to be like a pulse train with zero phase in the center as long as the input isn't shifted. For the sinc function, I think the phase ought to just be flat at zero for the width of the tophat in the magnitude spectrum but may have floating point junk outside. For the Gaussian, that seems to behave nicely and have flat phase of zero where the magnitude is high. Code below:

clc clearvars close all tSamp = 0.01; fSamp = 1/tSamp; t = -2:tSamp:2 - tSamp; % t = -2:tSamp:2; pick = 1; % 1 for rect, 2 for sinc, 3 for gaussian phaseCorrect = true; timeShift = 0; padFlag = 0; % 1 for end only, 2 for symmetric nSamp = length(t); % number*2=even % Want this to be able to handle signals with arbitrary delay, but not % getting expected results for zero-centered signals % x = rectpuls(t, 1); x = (abs((t-timeShift)) <= 0.5); y = 3*sinc(3*(t-timeShift)); sigma = 0.025; z = exp((-(t-timeShift).^2)/sigma^2); sigs = [x; y; z]; figure plot(t, sigs(pick, :)) xlabel('Time (s)') ylabel('Amplitude') if padFlag == 1 nSamp = 2*length(t); % number*2=even padSize = nSamp - length(t); pad = [sigs(pick, :), zeros(1, padSize)]; elseif padFlag == 2 nSamp = 2*length(t); % number*2=even padSize = nSamp - length(t); pad = [zeros(1, padSize/2), sigs(pick, :), zeros(1, padSize/2)]; else nSamp = length(t); pad = sigs(pick, :); end discSpec = fft((pad)); freq = (0:nSamp-1)/tSamp/nSamp; figure plot(abs(discSpec)) xlabel('Index') ylabel('Amplitude') title('Uncorrected Spectrum') % Dividing by fSamp=multiplying by tSamp to go from summation to % (approximate) integration % Multiplying by complex exponential to correct for the fact that the DFT % assumes the first sample is at t=0 and these start at t=-2. freq2 = -fSamp/2:fSamp/(nSamp):fSamp/2 - fSamp/(nSamp); tShift = t(1); shifter = exp(-j*2*pi*tShift*freq2); if phaseCorrect approxSpec = shifter.*fftshift(discSpec)/fSamp; else approxSpec = fftshift(discSpec)/fSamp; end figure freq2 = -1/(2*tSamp):1/(nSamp*tSamp):1/(2*tSamp) - 1/(nSamp*tSamp); plot(freq2, (abs(approxSpec))) xlabel('Frequency (Hz)') ylabel('Amplitude') title('Scaled Amplitude Spectrum') figure plot(freq2, (unwrap(angle(approxSpec)))) xlabel('Frequency (Hz)') ylabel('Phase') title('Unwrapped Phase Spectrum') figure plot(freq2, ((angle(approxSpec)))) % plot(freq2, atan2(imag(approxSpec), real(approxSpec))) xlabel('Frequency (Hz)') ylabel('Phase') title('Phase Spectrum') figure % plot(freq2, unwrap((angle(exp(j*2*pi*2*freq2))))) plot(freq2, unwrap(atan2(imag(shifter), real(shifter))))

EDIT 1: The primary test case I am trying to replicate is to get an accurate spectrum for $x(t)=2*sinc(3.2(t-1.5))$ where the sample spacing is 0.01 and $ -1\leq t \leq 4 $.
I think I follow the interpolation idea for a case where there are no restrictions placed, but perhaps my expectations for the result are out of order. For this case, the phase spectrum should have a negative slope of about $2\pi(1.5)$ and be zero outside the region, right?

Improve this question
$\endgroup$
5
  • $\begingroup$ Welcome to dsp.stackexchange. I love your name, hopefully we can change that and see you back here as "enlightened" $\endgroup$ Commented Oct 11 at 12:58
  • $\begingroup$ Hi Stuck- I see your edit which I didn't address in my answer. And I suspect your assignment is more about understanding phase in frequency versus delay in time which now explains the step 2 in your question above. So my answer makes the magnitude spectrum match but DOES NOT ADDRESS THE PHASE. How I would proceed so as to not cause and error (and allow you to do the work) is start with a zero phase condtion (2 * sinc(3.2t)) symetric about the time axis so extend -2.5 to 2.5. Do what you need per your step 2 to ensure it is 0 phase. Then carefully shift time in your function and in the range $\endgroup$ Commented Oct 11 at 16:18
  • $\begingroup$ And in doing that review the fundamental Fourier properties that relate a shift in time to a phase slope in freq (your step 2) to make sure you are moving in the right direction and to the right scale. This is a worth while exercise for you to work through slowly and carefully, just let us know if you're still "Stuck". ( I suspect you could complete your assignment by addressing the phase alone with the original samples you have and commenting that the magnitude is off due to frequency domain aliasing and what would need to be done to address that- my suspicion being this was all about phase) $\endgroup$ Commented Oct 11 at 16:20
  • $\begingroup$ @DanBoschen I think my primary issue is that I'm not sure how to interpret the phase plot, so I will try making another question to focus in on that. It does, however, sound like you agree that the correction methods I'm using are accurate, so is that the case? $\endgroup$ Commented Oct 11 at 16:50
  • $\begingroup$ Yes and good plan to make another question specific to phase, and see if you understand the FT of $x(t-\tau)$ is $e^{-j\omega \tau}$ and how to interpret that. $\endgroup$ Commented Oct 12 at 0:42

2 Answers 2

Reset to default
0
$\begingroup$

The difference between the Continuous Time Fourier Transform and the Discrete Fourier Transform (the DFT, which the FFT is an algorithm to compute) is frequency domain aliasing (assuming the time domain waveform as truncated to the fixed time duration used is the Fourier Transform we seek, otherwise truncation errors are another consideration and brings in the details of windowing).

To estimate the true Fourier Transform using the DFT, simply resample or interpolate the time domain waveform to a higher rate until the frequency domain aliasing is insignificant. This means, as in A/D conversion, meet Nyquist's criterion by sampling at higher than twice the occupied bandwidth of the signal.

Thus we see the impossibility with the concept of a rectangular pulse in either time or frequency: The Fourier Transform of a rectangular pulse is a Sinc function, so if we had a rectangular pulse in time, the resulting Sinc function in frequency extends (with relatively slow decay) to positive and negative infinity. In order to possibly meet the Nyquist criterion stated, it is necessary for that pulse to transition slowly (eliminating the high frequency components). Once we do that, then the FFT of the modified pulse in time will be nearly identical to samples of the continuous time Fourier Transform in frequency. (If you want to see more samples in frequency of that same CTFT, then simply zero pad the waveform as well prior to taking the FFT).

This is very similar to the process I use to make high quality filters, which requires approximating the true inverse Fourier Transform from the inverse DFT (and then windowing the time domain function after that). The first part prior to the windowing is similar to what the OP is asking for here in reverse (and applies given the reciprocity of the Fourier Transform). The approach I outline is very simple and convenient when processing resources are not an issue (so post processing applications rather than real-time hardware implementations where we may care otherwise). For further details including how to make great FIR filters once the true time domain impulse response is known (inverse Fourier Transform), please see this recent blog post:

https://dsprelated.com/showarticle/1760.php

Improve this answer
$\endgroup$
0
$\begingroup$

Adding to Dan's response.

Figure 2.2 from Digital Spectral Analysis by S. Lawrence Marple is a handy reference for relating the Continuous and Discrete Time Transforms. They’re related by windowing and sampling operations in both domains. enter image description here

The windowing and sampling functions of one domain have an equivalent operation in the other. TS is equivalent to periodic extension in frequency (frequency replication), FS periodic extension in time. TW is equivalent to convolution with sinc() in frequency, and FW is sinc() convolution in time.

Even though both are computed from discrete time data, notice the distinction between the Discrete Time Fourier Transform (DTFT) which is continuous in frequency and the Discrete Time Fourier Series (DTFS) which is calculated at discrete frequencies. The FFT, an algorithm for computing the DTFS, operates on a finite length record of discrete time samples and computes outputs at discrete frequencies.

Understanding these relationships will help you write a better CTFT approximation using DTFS.

Definitions:

CTFT $\begin{equation}X(f)=\int_{-\infty }^{\infty}x(t)e^{-j2\pi ft}dt\end{equation}$

DTFT $X(f) = T\sum_{n=-\infty}^{\infty}x[n]e^{-j2\pi fnT}$ $-1/{2T} \le f\le 1/{2T}$

DTFS $X[k] = T\sum_{n=0}^{N-1}x[n] e^{-j2\pi kn/N}$ $k=-\frac{N}{2}, ...\frac{N}{2}-1$

CTFS $X[k]=\int_{0 }^{NT}x(t)e^{-j2\pi kt/NT}dt$ $k=0, \pm 1, \pm 2, ..., \pm \infty$

Improve this answer
$\endgroup$

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.