View Full Version : DFT in C
Hello to everyone :) Since I did not notice a presentation section I put some info about myself before the question,.
My name is Enrico, I am Italian (form Modena) and I am 18. I am studying Electronics&Telecommunications in High School and I have applied for UCSD/UCLA/Berkeley (I hope they want me :smile: ). If you have questions, feel free to ask more about me.
Here the question. I am working on a routine which should perform the DFT of a discrete signal in C. More or less it works, but I noticed that the amplitude of the harmonics is not right. Here a little plot of the spectrum (is this the right word?) af a square wave, whose amplitude is (should be) 255.
http://img407.imageshack.us/img407/4979/quadramigliore.jpg
Here there is my program:
//---------------------------------------------------------------------------
#include <clx.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#define N 2048
#pragma hdrstop
//---------------------------------------------------------------------------
#pragma argsused
int main(int argc, char* argv[])
{
FILE *fp;
int Input[N];
int OutputR[N/2];
int OutputI[N/2];
int OutputAmp[N/2];
int OutputPha[N/2];
double Hamming[N];
int i,k;
double angolo;
//input square wave
for(i=0;i<N;i++){
angolo = (i*M_PI)/180;
Input[i] = (int)(324*sin(angolo) + 108*sin(3*angolo) + 65*sin(5*angolo) + 47*sin(7*angolo) + 36*sin(9*angolo) + 30*sin(11*angolo) + 25*sin(13*angolo) );
}
fp = fopen("file.txt","w");
for(i=0;i<N;i++){
fprintf(fp,"%d \n",Input[i]);
}
//Hamming Window data
for(i=0;i<N;i++){
Hamming[i] = 0.54-0.56*cos(i*M_PI/512);
Input[i] = (int)((double)Input[i]*Hamming[i]);
}
//DFT
//azzero i vettori
for(i=0;i<N/2;i++){
OutputR[i] = 0;
OutputI[i] = 0;
}
//calcolo i valori con la formula
for(k=0;k<N/2;k++){
for(i=0;i<N;i++){
OutputR[k] = OutputR[k]+Input[i]*cos(2*M_PI*k*i/N);
OutputI[k] = OutputI[k]-Input[i]*sin(2*M_PI*k*i/N);
}
}
//conversione polare->rettangolare
for(i=0;i<N/2;i++){
OutputAmp[i] = (int)sqrt(abs(OutputR[i]^2+OutputI[i]^2));
}
//scrivo su file
fp = fopen("file1.txt","w");
for(i=0;i<N/2;i++){
fprintf(fp,"%d \n",OutputAmp[i]);
}
fclose(fp);
return 0;
}
//---------------------------------------------------------------------------
Does someone have advice?
Thank you.
Hello to everyone :) Since I did not notice a presentation section I put some info about myself before the question,.
My name is Enrico, I am Italian (form Modena) and I am 18. I am studying Electronics&Telecommunications in High School and I have applied for UCSD/UCLA/Berkeley (I hope they want me :smile: ). If you have questions, feel free to ask more about me.
Here the question. I am working on a routine which should perform the DFT of a discrete signal in C. More or less it works, but I noticed that the amplitude of the harmonics is not right. Here a little plot of the spectrum (is this the right word?) af a square wave, whose amplitude is (should be) 255.
http://img407.imageshack.us/img407/4979/quadramigliore.jpg
Here there is my program:
I added [ code] and [ /code] tags around your code so that your indentation will be preserved.
//---------------------------------------------------------------------------
#include <clx.h>
#include <stdio.h>
#include <conio.h>
#include <math.h>
#define N 2048
#pragma hdrstop
//---------------------------------------------------------------------------
#pragma argsused
int main(int argc, char* argv[])
{
FILE *fp;
int Input[N];
int OutputR[N/2];
int OutputI[N/2];
int OutputAmp[N/2];
int OutputPha[N/2];
double Hamming[N];
int i,k;
double angolo;
//input square wave
for(i=0;i<N;i++){
angolo = (i*M_PI)/180;
Input[i] = (int)(324*sin(angolo) + 108*sin(3*angolo) + 65*sin(5*angolo) + 47*sin(7*angolo) + 36*sin(9*angolo) + 30*sin(11*angolo) + 25*sin(13*angolo) );
}
fp = fopen("file.txt","w");
for(i=0;i<N;i++){
fprintf(fp,"%d \n",Input[i]);
}
//Hamming Window data
for(i=0;i<N;i++){
Hamming[i] = 0.54-0.56*cos(i*M_PI/512);
Input[i] = (int)((double)Input[i]*Hamming[i]);
}
//DFT
//azzero i vettori
for(i=0;i<N/2;i++){
OutputR[i] = 0;
OutputI[i] = 0;
}
//calcolo i valori con la formula
for(k=0;k<N/2;k++){
for(i=0;i<N;i++){
OutputR[k] = OutputR[k]+Input[i]*cos(2*M_PI*k*i/N);
OutputI[k] = OutputI[k]-Input[i]*sin(2*M_PI*k*i/N);
}
}
//conversione polare->rettangolare
for(i=0;i<N/2;i++){
OutputAmp[i] = (int)sqrt(abs(OutputR[i]^2+OutputI[i]^2));
}
//scrivo su file
fp = fopen("file1.txt","w");
for(i=0;i<N/2;i++){
fprintf(fp,"%d \n",OutputAmp[i]);
}
fclose(fp);
return 0;
}
//---------------------------------------------------------------------------
Does someone have advice?
Thank you.
Is there some reason all of your arrays except the Hamming array are type int? You might be losing some precision by not making them double.
Also, you should double-check the coefficients that you are using in your Fourier calculation.
Thenk you for the answer :)
I chose to use only int to "simulate" the output of an ADC. I would like to write a FFT/DFT software for a Z80 processor, but beore I wanted to see what happened with integer value on the DFT. Maybe is better the output arrays to be double, isn't it? I make the change right now.
I also checked again the Fourier coefficients and they are right. I use this formula:
Y_{n} = \frac{4\cdot Y_{M}}{n\cdot \pi }
n = 1,3,5, ...
EDIT: I can not use double in Output arrays because the compiler says "Illegal use of floating point" in the conversion polar to rectangular.
I assume you mean this code.
//conversione polare->rettangolare
for(i=0;i<N/2;i++){
OutputAmp[i] = (int)sqrt(abs(OutputR[i]^2+OutputI[i]^2));
}
Instead of the abs() function in math.h, which takes an int or long value and returns an int or long, use the abs() function declared in stdlib.h. See http://www.cplusplus.com/reference/clibrary/cstdlib/.
Yes I meant that part of code.
How can I choose which abs() function to use? They have the same name(in math.h and stdlib.h) ... Other advice? I think that maybe the Hamming window is not totally right, but I do not know how to make it better.
The compiler will know which version of abs() to use based on the arguments. Since math.h and stdlib.h include versions of abs() with different signatures (i.e., different parameter types), if you call abs() with a double argument, the compiler will pick the overloaded version with that parameter type.
Can you be more specific about why you think your Hamming window code isn't right?
The compiler will know which version of abs() to use based on the arguments. Since math.h and stdlib.h include versions of abs() with different signatures (i.e., different parameter types), if you call abs() with a double argument, the compiler will pick the overloaded version with that parameter type.
Can you be more specific about why you think your Hamming window code isn't right?
Ok I will try to use the right abs().
I am not completely sure of the constants used in the Hamming Window.
Hamming[i] = 0.54-0.46*cos(i*M_PI/512);
0.54 and 0.46 are OK. I found the formula on the web and taking a look at the plot of the function 512 seemed a reasonable value for M. How can I find the best M ?
What about the amplitudes of the harmonics, which are not right?
Sorry, I don't know anything about the Hamming window, so I can't suggest a reasonable value for M.
The compiler will know which version of abs() to use based on the arguments. Since math.h and stdlib.h include versions of abs() with different signatures (i.e., different parameter types), if you call abs() with a double argument, the compiler will pick the overloaded version with that parameter type.
Can you be more specific about why you think your Hamming window code isn't right?
As of C99, there is no overloading in C. C1x has plans for a preprocessor based method for making type generic expressions, which can be used to spoof overloading. You're thinking C++.
Hamming[i] = 0.54-0.46*cos(i*M_PI/512);
What you should have is
Hamming[i] = 0.54 - 0.46 * cos(i*M_PI/N);
The window width needs to match the DFT width.
Also in math.h is the function double fabs(double);
I had time to take a look at my routine. First, the problem was not only in the abs function (by the way fabs() works well, thank you DrGreg ;) ), but seems that i cannot square a double in this way variable^2, I have to use the function pow(x,y). Using this function and doubles instead of int does not work, whe I plt the graph the number are too high (the smaller one is around 100.000.000). I also tried to match the window width with the width of the data, but the spectrum is worst.
Also in math.h is the function double fabs(double);
I was thinking that there was an fabs function that should be used instead of abs, but found a Web page that talked about an abs overload.
I had time to take a look at my routine. First, the problem was not only in the abs function (by the way fabs() works well, thank you DrGreg ;) ), but seems that i cannot square a double in this way variable^2, I have to use the function pow(x,y). Using this function and doubles instead of int does not work, whe I plt the graph the number are too high (the smaller one is around 100.000.000). I also tried to match the window width with the width of the data, but the spectrum is worst.
I've gotten so used to working with LaTeX that I did spot those variable^2 expressions. There is a ^ operator in C/C++, but it is the exclusive or operator, not an exponent operator. This is guaranteed to give you bad results. A better choice than using the power function is to just use multiplication, like so:
OutputR[i] * OutputR[i] + OutputI[i] * OutputI[i])
Thank you so much :) Now the plot of the spectrum is really clear ... just beautiful!
But the problem still remains ...
http://img713.imageshack.us/img713/3172/spettroscalasbagliata.jpg
If you take a look at the graph you can see that it is a "good" spectrum qualitative, but a wrong quantitative one, because the amplitudes of the harmonics are all wrong ...
Does anyone know hot so solve the problem?
vBulletin® v3.8.7, Copyright ©2000-2012, vBulletin Solutions, Inc.