0
$\begingroup$

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 ?

Improve this question
$\endgroup$
3
  • $\begingroup$ 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". $\endgroup$
    – Cedron Dawg
    Commented Jul 19, 2019 at 14:31
  • $\begingroup$ @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 ? $\endgroup$
    – SheppLogan
    Commented Jul 19, 2019 at 15:30
  • $\begingroup$ 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. $\endgroup$
    – Cedron Dawg
    Commented Jul 19, 2019 at 15:46

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
17
  • $\begingroup$ 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? $\endgroup$
    – SheppLogan
    Commented Jul 19, 2019 at 11:49
  • $\begingroup$ 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? $\endgroup$
    – SheppLogan
    Commented Jul 19, 2019 at 11:54
  • $\begingroup$ @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. $\endgroup$
    – Florian
    Commented Jul 19, 2019 at 12:24
  • $\begingroup$ 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) ? ;) $\endgroup$
    – SheppLogan
    Commented Jul 19, 2019 at 12:35
  • 1
    $\begingroup$ In regards of DFT, the dot-product and correlation I really enjoyed the following tutorial: SEEING CIRCLES, SINES, AND SIGNALS $\endgroup$
    – Irreducible
    Commented Jul 19, 2019 at 12:42

Your Answer

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

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