3
$\begingroup$

I got an $N$, in my case 512, point FFT of a real-valued signal. Based on some calculation in my application I determine the parameters $k \in [1, N-1]$, the number of oscillations per period, $\phi \in [0, 2\pi[$, phase. I use these parameters and the function $\cos\left(2\pi\cdot \frac{k}{N}\cdot t + \phi\right)$ to describe a signal in time-domain. Instead of calculating all the points of the signal and then doing an FFT I would like to directly calculate the result in frequency domain.

I calculated the FFT coefficients for the function as following:

$$c_n = \frac{N}{2}\left(\mathrm{e}^{i\phi} \cdot \mathrm{sinc}\left(n-k\right)+\mathrm{e}^{-i\phi} \cdot \mathrm{sinc}\left(n+k\right)\right), \quad\text{ with }\mathrm{sinc}(x) = \frac{\sin(\pi x)}{\pi x}$$

Note: $N$ is in here because the output of the FFT algorithm I use is not scaled and scaling is done in IFFT:

I verified my coefficients using Mathematica with the FourierCoefficient function.

I calculated a 512 elements long vector of complex doubles in Matlab and set my first vector to zero as there is no (dc) offset.

  • Is this valid or do I have to calculate bin 0 too if my $k$ is not a whole number?
  • Can I use $c_0$ for it?

The next 255 values I filled $c_n$ with $n$ from -1 to -255. The next value I filled with zero as I was unsure which coefficient should go there. And then the next 255 values I filled with $c_n$ with $n$ from 255 to 1.

To verify the correctness of my generated signal I plotted the real value of the IFFT of my vector. The output was as expected as long as $\phi$ was 0 and $k$ was a whole number.

Plot (blue is what was expected and red the result from IFFT) with $\phi=0$ and $k = 8.3$: Plot with phi=0 and k=8.3

Plot with $\phi=0.3$ and $k = 8$: Plot with phi=0.3 and k=8

Plot with $\phi=0.3$ and $k = 8.3$: Plot with phi=0.3 and k=8

For reference my matlab code is:

k = 8.3; %number of oscillations per period
phi = 0.3; %phase in radians

t = linspace(0,511,512);
n = linspace(-1,-255,255);
fft_coeff0 = 256*(exp(i*phi)*sinc(-k)+exp(-i*phi)*sinc(k));
fft_coeffs(1:255) = 256*(exp(i*phi)*sinc(n-k)+exp(-i*phi)*sinc(n+k));

reconst_fft = [fft_coeff0 fft_coeffs 0 conj(fft_coeffs(end:-1:1))];

y = cos(k/512*2*pi*t+phi);
plot(t, ifft(reconst_fft), 'r', t, y, 'b')
  • What could cause the difference between my expected output and the IFFT?
  • How can I fix it?

EDIT: Changed code to only calculate one side and use complex conjugating for the other side and added $c_0$ coefficient to vector.

Plot with $\phi=0.0$ and $k = 8.5$: Plot with phi=0.3 and k = 8.3:

Adding $c_0$ fixed the wrong amplitude.

$\endgroup$
3
  • $\begingroup$ I am not sure if the sinc function is correct in your calculation of $c_n$. I would expect something like $\frac {\sin(\pi (k-n))}{\sin(\pi (k-n)/N)}\exp(j\pi (k-n)(N-1)/N)$ instead of $sinc(\pi(k-n))$. This should come from the DFS of the windowing. $\endgroup$
    – Harris
    Sep 1, 2015 at 18:13
  • $\begingroup$ @Harris You're right. I derived my $c_n$ from the continuous fourier series. Nevertheless if using the DFS for calculating for large $N$ the error will vanish. I modified my matlab script to test with the DFS formula which slightly changed the output but didn't solve my actual problem. $\endgroup$
    – wiizzard
    Sep 2, 2015 at 22:42
  • $\begingroup$ Another problem is the fractional carrier frequency. If this is not a rational number then your signal is not really periodic. $\endgroup$
    – Harris
    Sep 3, 2015 at 16:57

3 Answers 3

2
$\begingroup$

If you have a pure tone signal defined by:

$$ S[n] = M \cos( \alpha n + \phi ) $$

Where $f$ is the frequency in cycles per frame and $\alpha$ is in radians per sample:

$$ \alpha = f \cdot \frac{ 2\pi }{ N } $$

And use a $ 1/N $ normalized DFT:

$$ Z[k] = \frac{ 1 }{ N } \sum_{n=0}^{N-1} { S[n] e^{ -i\beta_k n } } $$

Where:

$$ \beta_k = k \cdot \frac{ 2\pi }{ N } $$

The most concise formula for calculating the DFT bin values is:

$$ Z[k] = \frac{ M }{ 2N } \left[ \frac{ Ue^{ i\beta_k } - V }{ \cos( \alpha ) - \cos( \beta_k ) } \right] $$

Where:

$$ U = \cos( \alpha N + \phi ) - \cos(\phi) $$

$$ V = \cos( \alpha N -\alpha + \phi ) - \cos( -\alpha + \phi ) $$

Notice that there is a pluggable discontinuity at the bin value of a whole integer frequency. When the frequency is very close to an integer value there is a more computationally robust form of the formula:

$$ Z[k] = \frac{ M }{ 2 } \cdot \frac{ \sin \left( \frac{ \delta N }{2} \right) }{ N \sin \left( \frac{ \delta }{2} \right) } \cdot \left[ \frac{ \sin \left( \frac{ \delta N }{2} + \phi \right) e^{ i\beta_k } - \sin \left( \frac{ \delta N }{2} - \delta - \beta_k + \phi \right) }{ \sin \left( \frac{ \delta }{2} + \beta_k \right) } \right] $$

Where:

$$ \delta = \alpha - \beta_k = ( f - k ) \cdot \frac{ 2\pi }{ N } $$

The derivations of these formulas can be found in my blog articles "DFT Bin Value Formulas for Pure Real Tones" and "An Alternative Form of the Pure Real Tone DFT Bin Value Formula".

I realize this is a reply to a very old post, but Community dredged it up and finding these formulas, and the subsequent inverses which give exact frequency formulas is what got me into DSP.

Ced

$\endgroup$
0
$\begingroup$

If you want to reconstruct a strictly real result sine wave, the input to the IFFT needs to be conjugate symmetric. Whereas your code appears to have not inverted the imaginary component in the upper/negative half of the FFT coefficients. C0 should from your equation, not zero, as non-periodic-in-window sinewaves can end up with a DC offset in that window.

$\endgroup$
1
  • 1
    $\begingroup$ I changed my code to only calculate one half but as Cn is symmetric it should not matter though. Adding C0 fixed the problem with the incorrect amplitude. $\endgroup$
    – wiizzard
    Apr 8, 2014 at 18:10
0
$\begingroup$

If I understand the question correctly, this result may be of help. Keeping in mind that amplitude modulation is convolution in frequency domain. All you need to do then is convolve the unit impulse in frequency domain by the step function:

$$H({\omega}) = e^{j \omega\frac{N-1}{2}}\frac{\sin\omega\frac{N}{2}}{\sin\frac{\omega}{2}}$$

The equation for the cosine wave is:

$$f[n] = \cos(\omega_0 n) \overset{DTFT}{\leftrightarrow} F(\omega)=1, \omega =\omega_0$$

And likewise for the negative side of the spectrum.

The convolution is just a frequency shifted version of the window: $$F(\omega) * H(\omega) = \frac{1}{2}\left (e^{j (\omega-\omega_0)\frac{N-1}{2}}\frac{\sin\{(\omega-\omega_0)\frac{N}{2}\}}{\sin\frac{(\omega-\omega_0)}{2}} + e^{j (\omega+\omega_0)\frac{N-1}{2}}\frac{\sin\{(\omega+\omega_0)\frac{N}{2}\}}{\sin\frac{(\omega+\omega_0)}{2}} \right) $$

Converting this to the DFT representation should not be too difficult. As was hinted in the comments, the rectangular window is not a sinc function in discrete time, I think that's one of your errors. Also I think you have one extra minus sign on the 2nd exponential.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.