0

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 ?

SheppLogan
  • 683
  • 8
  • 22
  • For a better understanding of how the tone parameters result in bin values, check out https://www.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 : https://dsp.stackexchange.com/questions/59305/using-fourier-coefficients-to-reconstruct-data-in-matlab/59306#59306

    X[k] and X[n-k] are the positive and negative freq bins for freq k, yes ?

     – SheppLogan 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: https://dsp.stackexchange.com/questions/59068/how-to-get-fourier-coefficients-to-draw-any-shape-using-dft $$ $$ 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

1 Answers1

2

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.

Florian
  • 2,463
  • 1
  • 6
  • 16
  • 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? – SheppLogan 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? – SheppLogan 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) ? ;) – SheppLogan Jul 19 '19 at 12:35
  • 1
    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
  • This tutorial look really cool, thank you very much - i will definitely read it this weekend, thanks. – SheppLogan Jul 19 '19 at 12:48
  • My background is signal processing for multichannel systems, in particular RF (radar/comms). I'm sure some people will insist correlation only applies to random processes. Yet, it's not uncommon to "abuse" the term for other things. It's a matter of convention which varies widely between disciplines. – Florian Jul 19 '19 at 13:04
  • nice,thanks so much. If you have time maybe you can clarify my "phase" problem in question https://dsp.stackexchange.com/questions/59594/dct-vs-dft-why-do-we-need-want-phase/59596?noredirect=1#comment119098_59596 ->last comment, why is phase = -pi/2 and pi/2 (for the negative resp. positive spikes ) for a pure sine single freq (e.g. sin(2pif*t) with f=2 in my example) – SheppLogan Jul 19 '19 at 14:15
  • @Machupicchu The term "correlation" is a slippery one around here . (echoing Florian) It has a different definition in DSP vs statistics. I would avoid the term, but embrace the concept. – Cedron Dawg Jul 19 '19 at 14:26
  • Ok, so I can use the term dot products and similarity, instead, I guess – SheppLogan Jul 19 '19 at 14:27
  • @Machupicchu, Yeah, you can do that. This answer may help you with an interpretation of what the bin values mean: https://dsp.stackexchange.com/questions/59305/using-fourier-coefficients-to-reconstruct-data-in-matlab/59306#59306 – Cedron Dawg Jul 19 '19 at 14:29
  • @CedronDawg thanks. The formula you give: $e^{i n (N - k) \frac{2\pi}{N} } = e^{-i n k \frac{2\pi}{N} }$ does it come from a property of complex exponentials? – SheppLogan Jul 19 '19 at 14:45
  • @Machupicchu Yes, very much so.

    $$ e^{i n (N - k) \frac{2\pi}{N} } = e^{i N k \frac{2\pi}{N} } \cdot e^{-i n k \frac{2\pi}{N} } = e^{i k 2\pi } \cdot e^{-i n k \frac{2\pi}{N} } = (e^{i 2\pi })^k \cdot e^{-i n k \frac{2\pi}{N} } = 1^k \cdot e^{-i n k \frac{2\pi}{N} }$$

    When $k$ is an integer. It is the reason for the periodicity of the DFT. See also: https://www.dsprelated.com/showarticle/754.php

     – Cedron Dawg Jul 19 '19 at 14:57
  • Thanks. And the X[k] and X[N-k] are they corresponding to the positive, respectively negative DFT coefficients for the frequency k ? if I understand correctly? Although in Matlab you can't do X(N-k) rather X(N-k+1-1) because of 1-based indexing I guess? – SheppLogan Jul 19 '19 at 15:06
  • @Machupicchu In MATLAB, which I don't use, the bin values for $k$ will be at X[k+1] and X[N-k+1]. We are getting the chat prompt, please put your next comment under your question, and we won't bother Florian either ;-). – Cedron Dawg Jul 19 '19 at 15:20
  • @Florian This time another question for you, - hope i m not bothering too much - I noticed that when plotting the phase i.e. plot(angle(X)) from your code above I get sort of a 'random' zigzag. I guess it is related to the rounding errors because phases for other than the +freq -freq bin since in the example its a pure sine wave - should not have anything to do, Therefore are they just rounding artifacts and we only care about the phases for X(3) and X(N-freq+1) ? – SheppLogan Jul 19 '19 at 15:40
  • For the number zero, its phase is undefined. When a number is "essentially zero", i.e., zero except rounding errors, its phase is pretty much random. So when plotting the phase, discard those values that have zero magnitude, the phase for them is not meaningful. In other words: yes, you may ignore them. – Florian Jul 19 '19 at 15:57