Calculation of cosine in frequency domain instead of calculatin in time-domain followed by a FFT

Question

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.3$ and $k = 8$:

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

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$:

Adding $c_0$ fixed the wrong amplitude.

Answer 1

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

Answer 2

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.

Answer 3

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.