Practical question about DFT in Matlab

Please consider this piece of code:

clc;clear all;
fs=10^3;
t=linspace(0,1,fs);
f=2;
x=sin(2*pi*f*t);

W=dftmtx(fs);

X=fft(x);

X(3)
X(end-1)

%check : imaginary parts are 180 degrees dephased but right frequency
%plot(imag(W(end-1,:)));hold on; plot(imag(W(3,:)))

Xrec=(1/length(X))*(X(end-1)*W(3,:)+X(3)*W(end-1,:));

plot(x,'+g');hold on;plot(Xrec,'r','LineWidth',1.25);


The idea is to show (for my personal understanding) that we can reconstruct the pure sine wave x (of frequency 2 cycles) with the only 2 Fourier coefficients (corresponding to the 2 spikes we see in the spectrum) i.e. a negative frequency = -2 and a positive frequency = 2, for this example.

It seems to work well with the line:

Xrec=(1/length(X))*(X(end-1)*W(3,:)+X(3)*W(end-1,:));


Where I want to show that the signal x can be reconstructed by the weighted contribution of W(3,:) and W(end-1,:) (which I found to correspond to +freq and -freq (+2/-2)). This is the point of view of linear algebra where the vector x is a linear combination of the basis vectors W(3,:) and W(end-1,:) with the complex weights X(3) and X(end-1).

Is this explanation correct?

Moreover, there is the perspective of "matched filter", i.e. correlation. From this point of view we should have only contribution of the 2 frequencies aforementioned here 2/-2. But I still get some non-zero correlation when doing e.g. :

W(555,:)*x' = -0.0063 - 0.0011i

... i.e. W(k != 3 or 999) (which are the indexes corresponding to freq=2 or freq=-2)

its small but not zero, so I would then to think that these have some contribution to the signal yet the should not ?

• For a better understanding of how the tone parameters result in bin values, check out dsprelated.com/showarticle/771.php. (19) gives the equation when you have a whole number of cycles in the frame. Otherwise, from (25) you can tell when $\alpha$ (the tone's radians per sample), is close to $\beta_k$ (Bin k's radians per sample) the magnitude of the bin gets larger. When your k's are further away, the magnitudes are smaller. This equation exactly describes "leakage". – Cedron Dawg Jul 19 '19 at 14:31
• @CedronDawg You mean append my next comments here ? Yes sorry of course X[N-k+1] the extra minus 1 didn't make sense So you in your other post : dsp.stackexchange.com/questions/59305/… X[k] and X[n-k] are the positive and negative freq bins for freq k, yes ? – Machupicchu Jul 19 '19 at 15:30
• Yes, a frequency of $-k$ and $N-k$ are aliases of each other. You can't tell them apart in the DFT. A striking example of what happens if you use the $N-k$ frequency instead of $k$ in your interpolations (Fourier series) can be seen visually for a complex signal with the "Fluffy Cloud" graphs in my answer here: dsp.stackexchange.com/questions/59068/…  If you don't interpret it as $-k$ you can't apply this $$\cos(\theta) = \frac{e^{i\theta}+e^{-i\theta}}{2}$$ for your Fourier series coefficient calculations. – Cedron Dawg Jul 19 '19 at 15:46

I augmented and "fixed" your code a little bit.

Your explanations are basically correct. Your main mistake is the following: your signal needs to be periodic in time domain, i.e., a whole number of periods must fit into your time window. If you use linspace this will not work: your first data point is 0, your last data point is 1, however cos(2*pi*f*0) = 1 = cos(2*pi*f*1). Hence your last and your first point are equal, they appear twice.

Try it with a very coarse spacing: cos(2*pi*linspace(0,1,5)) = [1, 0, -1, 0, 1]. Periodic copies give [1,0,-1,0,1,1,0,-1,0,1]. What you want is [1,0,-1,0,1,0,-1,0,1] though and for this you'd need to remove the last point.

A better way of doing it is this:

clear;
fs = 10^3;     % sampling frequency (Hz)
t0 = 1/fs;     % sample spacing in time (s)
NSamp = 1000;  % number of samples
f0 = fs/NSamp; % spacing in frequency
t_axis = (0:NSamp-1)*t0; % time axis
f_axis = (0:NSamp-1)*f0; % frequency axis
f_axis_symm = (-NSamp/2:NSamp/2-1)*f0; % frequency axis (symmetric)

f = 2;
x = sin(2*pi*f*t_axis);
X=fft(x);
C_two = X(2 + 1); % C[2]: need to add one since Matlab counts from one
C_minus_two = X(-2 + NSamp + 1); % C[-2]: need to add N (periodicity) and one (s.a.)

figure(1);
clf;
% stem(f_axis,abs(X)); % "asymmetric version [0,fs)
stem(f_axis_symm,fftshift(abs(X))); % symmetric" version [-fs/2,fs/2)
xlabel('f [Hz]');

W = conj(dftmtx(fs))/NSamp; % IDFT matrix: see help DFTMTX
Xrec = W(:,3)*C_two + W(:,NSamp-1)*C_minus_two; % reconstruct

figure(2);
clf;
plot(t_axis,x,'s','MarkerIndices',1:10:NSamp);
hold on;
plot(t_axis,Xrec,'b--');
xlabel('t [s]');


The code also shows the spectrum, plotted in a symmetric way so that you can see exactly two lines popping up at -2 Hz and 2 Hz. It also shows how to reconstruct the signal. Note that you need an IDFT matrix for reconstruction. See help dftmtx for a comment regarding this.

• Fantastic answer thanks. I had not thought about this periodicity wrt number of points. Now when I project (scalar product) any W(k != indexes of f= 2 or f=-2) I do get something of the order of 1e-17 which I guess I can simply consider =0 due to rounding error right? – Machupicchu Jul 19 '19 at 11:49
• An for the "theoretical" part, do you agree that I may say that the Fourier transform coefficients (in capital X) are the result of the "correlation" of lowercase x with a given frequency., i.e. that a high value of X represents a high "correlation" with this freq, (and a 0 or very low means not correlated to the other frequencies ),; In this example only two and minus_two are correlated, all others are 0 correlated - although the term "correlation" maybe is not the correct one here - what would you say? – Machupicchu Jul 19 '19 at 11:54
• @Machupicchu: Yes, errors of 1e-17 are just rounding errors. You can consider them zero. And yes, I think the DFT can be seen as correlating a given signal with harmonics of different frequencies. A large coefficient tells us that the signal is very similar to that given harmonic. Depends how you define correlation, some people associate it to randomness. I can assure you that it is also not uncommon to use the term in this context. If you want to avoid correlation: it's an inner product really. Large inner product = strong similarity. – Florian Jul 19 '19 at 12:24
• Great answer thanks. I don't really want to avoid it. But someone (with a phd!) once implied that wasn't the correct term to use (myself starting a phd now, I want to be as precise as possible with any terms I will use). What term would you use personally given your background (if I may ask) ? ;) – Machupicchu Jul 19 '19 at 12:35
• In regards of DFT, the dot-product and correlation I really enjoyed the following tutorial: SEEING CIRCLES, SINES, AND SIGNALS – Irreducible Jul 19 '19 at 12:42