DFT MathLab is wrong or I am wrong

[btb4198]

so this is from math lab

x is a 2 kHz sinewave

>> n = 0:1023;
>> t = n/8000;
>> x = sin(2000*2*pi*t);
>> X = fft(x);
>> plot(n,abs(X));

I have uploaded the picture of the graph ...
but is it wrong
2000 is the only one that should have a real value right ?

 

[ZetaGuy]

I think it is correct. The Fourier transform of sin(a t) is a sum of two Dirac deltas up and down from the base frequency a. Your plot shows the discretization of those deltas at 1/4 and 3/4 of the range.

I imagine this sounds confusing because the sine is a "pure tone" and the power spectrum should have only one peak. The power spectrum of a real signal is an even function, so in practice one only needs to show half of the spectrum, and it is that half that has one peak.

 

[jasonRF]

Looks reasonable to me. recall sin(x)=\left( e^{ix} - e^{-ix} \right)/2i, so a real sinusoid will have both a positive and negative frequency component. You should be able to compute which frequency bins correspond to +/- 2000 Hz, given your sample rate and the number of time samples you used.

jason

EDIT: looks like someone beat me to this!

 

[btb4198]

I think it is correct. The Fourier transform of sin(a t) is a sum of two Dirac deltas up and down from the base frequency a. Your plot shows the discretization of those deltas at 1/4 and 3/4 of the range.

I imagine this sounds confusing because the sine is a "pure tone" and the power spectrum should have only one peak. The power spectrum of a real signal is an even function, so in practice one only needs to show half of the spectrum, and it is that half that has one peak.

where do you get 1/4 and 3/4 from ?

and if sine is a "pure tone" ( like ringing a bell) why does it have two frequencies and how do you know the other frequency should be ?

 

[btb4198]

Looks reasonable to me. recall sin(x)=\left( e^{ix} - e^{-ix} \right)/2i, so a real sinusoid will have both a positive and negative frequency component. You should be able to compute which frequency bins correspond to +/- 2000 Hz, given your sample rate and the number of time samples you used.

jason

EDIT: looks like someone beat me to this!

so is that why there are two frequencies ?
one Positives and one Negative ?

so when someone does F(t) = sin(2000*2πt);
there are really two frequencies : one positive and one negative?
So is 2000 not really the real frequency ?
if so what are the real frequencies and how do you solve for them?

 

[jasonRF]

A question: the code you posted looks like matlab. Is mathlab different, or is it a typo?

Also, what reading have you done on discrete signals and systems in general, and DFT/FFT in particular? Your other DFT question a few days ago was essentially the same...

As I showed above, a real sinusoid of frequency 2000 Hz has a DFT with peaks at +2000 Hz and -2000 Hz. To map the bins of the FFT output to frequency depends upon the particular FFT software you are using (different packages pack the positive and negative frequencies differently). I do know that for matlab, the first N/2 are the positive, and the next N/2 are the negative. There is a command fftshift that makes 0 frequency (DC) the middle bin. So in matlab, for your example, you can plot
>> fs = 8000; %your sample rate
>> N = 1024; %the number of samples you are using
>> F = fs*[-N/2:1:N/2-1]/N;
>> plot(F, abs(fftshift(X)));

jason

 

[jasonRF]

so is that why there are two frequencies ?
one Positives and one Negative ?

so when someone does F(t) = sin(2000*2πt);
there are really two frequencies : one positive and one negative?
So is 2000 not really the real frequency ?
if so what are the real frequencies and how do you solve for them?

Yes, see my example above with sin(x). The DFT is essentially expanding your input signal in a series of complex sinusoids. So a signal that is exp(i w t) can be represented by a signle complex exponential, while sin(w t) is a combination of exp(i w t) and exp(-i w t), so needs two of the complex exponentials to represent it.

 

[btb4198]

A question: the code you posted looks like matlab. Is mathlab different, or is it a typo?

Also, what reading have you done on discrete signals and systems in general, and DFT/FFT in particular? Your other DFT question a few days ago was essentially the same...

As I showed above, a real sinusoid of frequency 2000 Hz has a DFT with peaks at +2000 Hz and -2000 Hz. To map the bins of the FFT output to frequency depends upon the particular FFT software you are using (different packages pack the positive and negative frequencies differently). I do know that for matlab, the first N/2 are the positive, and the next N/2 are the negative. There is a command fftshift that makes 0 frequency (DC) the middle bin. So in matlab, for your example, you can plot
>> fs = 8000; %your sample rate
>> N = 1024; %the number of samples you are using
>> F = fs*[-N/2:1:N/2-1]/N;
>> plot(F, abs(fftshift(X)));

jason

Yes that code is Matlab but I did not write it and I do not have matlab.
sorry back in school I would called it math lab...
I got it from the file I have attached.

I am doing the DFT by hand and with a C# program I am writing

 

[btb4198]

ok I ran that example with my code and I got
img F[255] = 108.542903447448
real F[255] = 109.044055059703
img F[256] = -325.630757652576
real F[256] = -325.131141522754

I just learned that to get frequencies you can to do

256/1023 = F/8000
F = 2001.95Hz that can is really close to 2000

I also did
(255/1023) =( F/8000)
F = 1994.13HZ
so I think it is save say go with the bigger number

also I did
F = (767/1023) * 8000 = 5998.04

Now how would I know that F = 5998.04Hz is wrong if I did know I pull 2000Hz in ?

 

[btb4198]

also,
right now
I have Bin Vs Frequency
how do I get amplitude vs frequency

 

[jasonRF]

ok I ran that example with my code and I got
img F[255] = 108.542903447448
real F[255] = 109.044055059703
img F[256] = -325.630757652576
real F[256] = -325.131141522754

I just learned that to get frequencies you can to do

256/1023 = F/8000
F = 2001.95Hz that can is really close to 2000

I also did
(255/1023) =( F/8000)
F = 1994.13HZ
so I think it is save say go with the bigger number

also I did
F = (767/1023) * 8000 = 5998.04

Now how would I know that F = 5998.04Hz is wrong if I did know I pull 2000Hz in ?

Note that you should be dividing by 1024 (the number of samples) instead of 1023. Also, do your indices in your code start with 0 or 1. That is, is the first frequency bin F[0] or F[1] ? If it is F[0] (which I suspect since C starts arrays at element 0) then the peak at F[256] means the frequency corresponds to
(256/1024)*8000 = 2000 Hz exactly, and the other peak is at
(768/1024)*8000 = 6000 Hz exactly.

In any case, I will let you sort out indexing and such.

Since you are sampling at 8000 Hz, your Nyquist is at 4000 Hz, hence your 6000 Hz really corresponds to 4000-6000 = -2000 Hz. This is the negative frequency part of your sinusoid. There is nothing wrong here, and it is unambiguous how to interpret your fft output. Again, what you are doing with your Fourier analysis is representing your initial time-domain signal as a finite sum of complex exponentials; the DFT is the way you compute the coefficients of that expansion (and the FFT is simply an algorithm that performs the DFT in a smart way to minimize computation). Since your signal is effectively sin(x) you expect to see both a exp(ix) and exp(-ix) term. If your input signal was simply exp(i*2*pi*2000*t) then you would only get a single peak instead of 2 peaks. You should try it and see!

jason

 

[btb4198]

Note that you should be dividing by 1024 (the number of samples) instead of 1023. Also, do your indices in your code start with 0 or 1. That is, is the first frequency bin F[0] or F[1] ? If it is F[0] (which I suspect since C starts arrays at element 0) then the peak at F[256] means the frequency corresponds to
(256/1024)*8000 = 2000 Hz exactly, and the other peak is at
(768/1024)*8000 = 6000 Hz exactly.

In any case, I will let you sort out indexing and such.

Since you are sampling at 8000 Hz, your Nyquist is at 4000 Hz, hence your 6000 Hz really corresponds to 4000-6000 = -2000 Hz. This is the negative frequency part of your sinusoid. There is nothing wrong here, and it is unambiguous how to interpret your fft output. Again, what you are doing with your Fourier analysis is representing your initial time-domain signal as a finite sum of complex exponentials; the DFT is the way you compute the coefficients of that expansion (and the FFT is simply an algorithm that performs the DFT in a smart way to minimize computation). Since your signal is effectively sin(x) you expect to see both a exp(ix) and exp(-ix) term. If your input signal was simply exp(i*2*pi*2000*t) then you would only get a single peak instead of 2 peaks. You should try it and see!

jason

you are right about the 0 index and I should use 1024 and not 1023
I too now got
2000
6000
how do you come up with N and Fs for any giving problem ?

is fs = nyquist theorem you sample twice the F of your incoming signal?

like if the incoming signal was a person's voice, the F would be between 80 Hz to 1100 Hz
so Fs would be 2200Hz right?

so would N be 1024?

how many samples to take, I do know the more you take the closer your anwer will be to the right F

 

[jasonRF]

You need to sample at at least twice the frequency of the highest signal. If you don't, the signals will end up in "aliased" frequency bins. The standard example is automobile wheels in movies - there are times when it looks like the wheel is going backwards, right? This is because of aliasing - either the frame-speed of the film isn't high enough, or the "frame speed" of the human visual system isn't fast enough. In your example the freuqency was 2000 and the sample rate 8000 so you are not aliased. You should try a sinusoid of 5000 Hz with your same sample rate, and see what happens.

jason

 

[ZetaGuy]

Assuming you have are using the same conventions as matlab does for the fast Fourier transform, then when Fourier transforming an array of real numbers, the result should satisfy x_k = conjugate(x_1024-k) for arrays that start at zero and have 1024 entries. The frequency calculation for the points above the middle mark should be as in (1024-768)/1024 * 8000 = 2000

If you are trying to read the fundamental frequency out of a noisy signal, there will be more to the task than just FFTs, and often no FFTs, as there are better algorithms for that.

 

[btb4198]

Assuming you have are using the same conventions as matlab does for the fast Fourier transform, then when Fourier transforming an array of real numbers, the result should satisfy x_k = conjugate(x_1024-k) for arrays that start at zero and have 1024 entries. The frequency calculation for the points above the middle mark should be as in (1024-768)/1024 * 8000 = 2000

If you are trying to read the fundamental frequency out of a noisy signal, there will be more to the task than just FFTs, and often no FFTs, as there are better algorithms for that.

like what?

 

[ZetaGuy]

How about some more details on the problem you want to solve? The literature on finding frequencies is vast and bit more information would help.

 

[btb4198]

I am reading in a signal from a Microphone
would that be too noisy for DFT?

 

[ZetaGuy]

In analyzing a microphone signal, I imagine you could be interested in understanding if it has any structure that can be explored by other algorithms, or you may be interested in seperating the sound into perceptual components, such as voice and a guitar, or you may be trying to reduce the noise in the signal.

Suppose you are intested in detecting the fundamental frequency of your speech, say to build an app that tells you whether you are in tune or not. In this case you would also have to deal with perceptual issues, such as pitch changing with intensity at higher frequencies, but any such algorithm would require as input the frequency of the signal.

Frequency is only defined for an infinitely long signal, such as sine wave. For finite signals one has to make an approximation. In doing a fast Fourier transform with 1024 samples, the hidden assumption is that those samples then repeat again and again. This hidden periodicity will introduce extreneous peaks in your transform. Window functions are used in conjunction with spectral analysis to limit the resulting spectrum to a region of interest and eliminate artifacts.

The signal from even a simple built in laptop microphone will provide enough data to measure voice pitch. The signal will be noisy, but it should be possible to work around it.

The algorithms for pitch detection could be organized into two broad categories: those that Fourier transform and use the frequency components and those that stay in the time domain. For algorithms in the frequency domain, search for Piszczalski who in the late 1970s developed a system of music transcription. There are today many improvements to his original algorithms. Other useful buzz terms to seach would be "cepstrum analysis" and "comb filters". In the time domain a basic algorithm is to count the average number of times the signal crosses its median value, or "zero crossing rate". Auto correlation methods and shifted curve matching are other common methods.

Many of the basic techniques are now in textbooks. The Wikipedia page on signal processing could be an entry point into the literature.