convert time domain info into frequency domain.

[mashhur]

Hello there,

I need an advise to convert time domain info into frequency domain (i need to know signal's strength in DB or mDB).

let me explain in a detail my issue.
i got 8192 (fifo size) digital numbers in 0.012 milliseconds and those values are VPP (voltage points).
so, in order to know signal'e strength i need to utilize FFT or Hilbert transform.

how to implement them in this situation.

i really appreciate if somebody help me out on this issue.
please post the url of C CODE if you know.

Thanks in advance!

 

[godfreyl]

If you only need to know the signal's strength in DB, then you don't need an FFT or Hilbert transform.

Can't you just import the data into a spreadsheet, calculate the average value, subtract that from all the datapoints to get the AC, then calculate the RMS value and convert to dB?

 

[thylacine1975]

Hiya mashhur

Your problem sounds like a perfect (N = a power of two, even!) problem to be solved with the FFTW ("Fastest Fourier Transform in the West" libraries, available at https://www.fftw.org/.

It's been a few years since I used their routines in anger, but I remember a) they were awesomely fast, and b) they were slightly convoluted to set up due to their multitude of possible configurations. I recall reading and re-reading the documentation a couple of times ... and then it worked as hoped! Use of the libraries first involves defining a "plan" descibing how the Fourier Transform is to be calculated, pointing it at input and output buffers and finally calling the fftw_execute(plan); function to let it do its work.

Although my coding is nothing to rave about about, let me include a standalone (console) application I wrote that called the libraries to give you a worked example to start from. This simple program was written to calculate power spectral densities from several-GB long raw binary files, sampled and written to disk by an Agilent 'Acqiris' digitiser.

Good luck!

---

Code C - [expand]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
/* acqSpect.c    JW 21 May 2009
 *
 *  Cuts an Acqiris data file into a number of segments, and performs a DFT on each 
 * segment.  The magnitude of the resulting complex data is then calculated to return the
 * power spectral density. The PSD's of each segment are accumulated (summed) to determine
 * the average distribution, which is written to disk.  (The DFT operation is performed
 * by a call to FFTW's ([url]www.fftw.org[/url]) double-precision DLL library.)
 *
 *  Usage:
 *    acqSpect [inputFile] [outputFile] [N]    - where N = Number of points of DFT
 *
 *  Compile with (for the PC platform): 
 *    mingw32-gcc -O2 -o acqSpect.exe acqSpect.c -L./ -llibfftw3-3 -lm
 *
 *  Revision history
 *  ----------------
 *
 *  v1.0 - 21 May 2009 - Created.
 *  v2.0 -  2 Jun 2009 - Modified to overlap (50%) and window (Welch) data segments.
 *                       Data that is beyond an integer multiple of segments is discarded.
 *  v2.1 - 26 Jun 2009 - Bug fix?  Switched to using fread and pointer casting to read
 *                       the int16 data from the Acqiris.  I'm not sure if I needed to, 
 *                       but there were some subtle hints of integer overflow/wraparound
 *                       weirdness in v2.0 that I never fully explored.
 */
 
#include <stdio.h>
#include "fftw3.h"
#include <stdlib.h>
 
int main(int argc, char *argv[])
   {
   unsigned char aqData[2];
   long i, N, totalSamples;
   int seg, nSegs;
   FILE *inFile, *outFile;
     double *in, *psd, *window;
   fftw_complex *out;
   fftw_plan plan;
   
   /* Parse input arguments */
   
   if (argc != 4) {
      printf("acqSpect v2.1 - JW 26/6/2009.\n");
      printf("Usage: acqSpect [inputFile] [outputFile] [N]\n\n");
      printf("\t inputFile = Acqiris data file (raw binary format)\n");
      printf("\t outputFile = ASCII file to write, length N/2 + 1 (Fourier frequencies 0..Fs/2)\n");
      printf("\t N = Number of points of DFT for calculation.\n");
      return(0);
      }
   
   sscanf(argv[3], "%ld", &N);
   N = 2*(N >> 1);
   if (N <= 0) {
      printf("N must be an even integer > 0.\n");
      return(1);
      }
 
   inFile = fopen(argv[1], "rb");  
   if (inFile == NULL) {
      printf("Cannot open Acqiris input file.\n");
      return(1);
      }
      
   outFile = fopen(argv[2], "w");  
   if (outFile == NULL) {
      printf("Cannot open output file for writing.  Out of disk space?.\n");
      fclose(inFile);
      return(1);
      }
 
   /* Examine input file */
   
   fseek(inFile, 0, SEEK_END);
   totalSamples = ftell(inFile)/2;
   fseek(inFile, 0, SEEK_SET);
   printf("%d input samples found.\n", totalSamples);
 
   nSegs = (2*totalSamples)/N - 1;
   if (nSegs == 0)  {
      printf("Input data length needs to be at least 1.5 x N...\n");
      return(1);
      }     
   
   /* Allocate space for DFT data and prepare FFTW's "forward plan" */
   
   printf("Pre-calculating FFT coefficients...\n");
   in = (double*) fftw_malloc(sizeof(double) * N);
   out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N/2+1));
   plan = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE); 
   psd = (double *)malloc(sizeof(double) * (N/2+1));
   window = (double *)malloc(sizeof(double) * N);
 
   if ((psd == NULL) || (in == NULL) || (out == NULL) || (window == NULL)) {
      printf("Ack!  Dynamic memory allocation failed!\n");
      return(1);
      } 
 
   for (i=0; i<(N/2+1); i++)   /* Zero output array */
      psd[i] = 0.0;
 
   for (i=0; i<N; i++)   /* Precalculate Welch window coefficients */
      window[i] = 1.0 - (2.0*(double)(i)/(N-1.0) - 1.0)*(2.0*(double)(i)/(N-1.0) - 1.0); 
 
   /* Start processing */
 
   printf("Calculation: 000 %% complete.");
   
   for (seg=0; seg<nSegs; seg++) {
      printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%3d %% complete.", 100*seg/nSegs + 1);
 
      fseek(inFile, seg*N, SEEK_SET);         /* Jump to start of segment */
      for (i=0; i<N; i++) {                   /* And convert LSB0 MSB0 LSB1 MSB1... integer data -> double precision */
         if (fread(aqData, 2, 1, inFile) != 1) {
            printf("\nRead of data failed from input file (at byte %ld).\n", ftell(inFile));
            return(1);
            } 
            
         in[i] = window[i]*(double)(*(short *)(aqData))/32768.0;
         }
 
      fftw_execute(plan);                     /* Calculate DFT */
      
      for (i=0; i<(N/2+1); i++)
         psd[i] += (out[i][0]*out[i][0]) + (out[i][1]*out[i][1]); /* Calculate abs()^2 */ 
      }
   
   /* Write output file */
   
   for (i=0; i<(N/2+1); i++)
      fprintf(outFile, "%E\n", psd[i]/(double)(nSegs));
 
   /* Housekeeping */
   
   printf("\nDone!\n"); 
   fclose(inFile);
   fclose(outFile);
 
   fftw_destroy_plan(plan);
   fftw_free(in); 
   fftw_free(out);
   free(psd);
   }

---

 

[mashhur]

Hiya mashhur

Your problem sounds like a perfect (N = a power of two, even!) problem to be solved with the FFTW ("Fastest Fourier Transform in the West" libraries, available at https://www.fftw.org/.

It's been a few years since I used their routines in anger, but I remember a) they were awesomely fast, and b) they were slightly convoluted to set up due to their multitude of possible configurations. I recall reading and re-reading the documentation a couple of times ... and then it worked as hoped! Use of the libraries first involves defining a "plan" descibing how the Fourier Transform is to be calculated, pointing it at input and output buffers and finally calling the fftw_execute(plan); function to let it do its work.

Although my coding is nothing to rave about about, let me include a standalone (console) application I wrote that called the libraries to give you a worked example to start from. This simple program was written to calculate power spectral densities from several-GB long raw binary files, sampled and written to disk by an Agilent 'Acqiris' digitiser.

Good luck!

---

Code C - [expand]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
/* acqSpect.c    JW 21 May 2009
 *
 *  Cuts an Acqiris data file into a number of segments, and performs a DFT on each 
 * segment.  The magnitude of the resulting complex data is then calculated to return the
 * power spectral density. The PSD's of each segment are accumulated (summed) to determine
 * the average distribution, which is written to disk.  (The DFT operation is performed
 * by a call to FFTW's ([url]www.fftw.org[/url]) double-precision DLL library.)
 *
 *  Usage:
 *    acqSpect [inputFile] [outputFile] [N]    - where N = Number of points of DFT
 *
 *  Compile with (for the PC platform): 
 *    mingw32-gcc -O2 -o acqSpect.exe acqSpect.c -L./ -llibfftw3-3 -lm
 *
 *  Revision history
 *  ----------------
 *
 *  v1.0 - 21 May 2009 - Created.
 *  v2.0 -  2 Jun 2009 - Modified to overlap (50%) and window (Welch) data segments.
 *                       Data that is beyond an integer multiple of segments is discarded.
 *  v2.1 - 26 Jun 2009 - Bug fix?  Switched to using fread and pointer casting to read
 *                       the int16 data from the Acqiris.  I'm not sure if I needed to, 
 *                       but there were some subtle hints of integer overflow/wraparound
 *                       weirdness in v2.0 that I never fully explored.
 */
 
#include <stdio.h>
#include "fftw3.h"
#include <stdlib.h>
 
int main(int argc, char *argv[])
   {
   unsigned char aqData[2];
   long i, N, totalSamples;
   int seg, nSegs;
   FILE *inFile, *outFile;
     double *in, *psd, *window;
   fftw_complex *out;
   fftw_plan plan;
   
   /* Parse input arguments */
   
   if (argc != 4) {
      printf("acqSpect v2.1 - JW 26/6/2009.\n");
      printf("Usage: acqSpect [inputFile] [outputFile] [N]\n\n");
      printf("\t inputFile = Acqiris data file (raw binary format)\n");
      printf("\t outputFile = ASCII file to write, length N/2 + 1 (Fourier frequencies 0..Fs/2)\n");
      printf("\t N = Number of points of DFT for calculation.\n");
      return(0);
      }
   
   sscanf(argv[3], "%ld", &N);
   N = 2*(N >> 1);
   if (N <= 0) {
      printf("N must be an even integer > 0.\n");
      return(1);
      }
 
   inFile = fopen(argv[1], "rb");  
   if (inFile == NULL) {
      printf("Cannot open Acqiris input file.\n");
      return(1);
      }
      
   outFile = fopen(argv[2], "w");  
   if (outFile == NULL) {
      printf("Cannot open output file for writing.  Out of disk space?.\n");
      fclose(inFile);
      return(1);
      }
 
   /* Examine input file */
   
   fseek(inFile, 0, SEEK_END);
   totalSamples = ftell(inFile)/2;
   fseek(inFile, 0, SEEK_SET);
   printf("%d input samples found.\n", totalSamples);
 
   nSegs = (2*totalSamples)/N - 1;
   if (nSegs == 0)  {
      printf("Input data length needs to be at least 1.5 x N...\n");
      return(1);
      }     
   
   /* Allocate space for DFT data and prepare FFTW's "forward plan" */
   
   printf("Pre-calculating FFT coefficients...\n");
   in = (double*) fftw_malloc(sizeof(double) * N);
   out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * (N/2+1));
   plan = fftw_plan_dft_r2c_1d(N, in, out, FFTW_MEASURE); 
   psd = (double *)malloc(sizeof(double) * (N/2+1));
   window = (double *)malloc(sizeof(double) * N);
 
   if ((psd == NULL) || (in == NULL) || (out == NULL) || (window == NULL)) {
      printf("Ack!  Dynamic memory allocation failed!\n");
      return(1);
      } 
 
   for (i=0; i<(N/2+1); i++)   /* Zero output array */
      psd[i] = 0.0;
 
   for (i=0; i<N; i++)   /* Precalculate Welch window coefficients */
      window[i] = 1.0 - (2.0*(double)(i)/(N-1.0) - 1.0)*(2.0*(double)(i)/(N-1.0) - 1.0); 
 
   /* Start processing */
 
   printf("Calculation: 000 %% complete.");
   
   for (seg=0; seg<nSegs; seg++) {
      printf("\b\b\b\b\b\b\b\b\b\b\b\b\b\b\b%3d %% complete.", 100*seg/nSegs + 1);
 
      fseek(inFile, seg*N, SEEK_SET);         /* Jump to start of segment */
      for (i=0; i<N; i++) {                   /* And convert LSB0 MSB0 LSB1 MSB1... integer data -> double precision */
         if (fread(aqData, 2, 1, inFile) != 1) {
            printf("\nRead of data failed from input file (at byte %ld).\n", ftell(inFile));
            return(1);
            } 
            
         in[i] = window[i]*(double)(*(short *)(aqData))/32768.0;
         }
 
      fftw_execute(plan);                     /* Calculate DFT */
      
      for (i=0; i<(N/2+1); i++)
         psd[i] += (out[i][0]*out[i][0]) + (out[i][1]*out[i][1]); /* Calculate abs()^2 */ 
      }
   
   /* Write output file */
   
   for (i=0; i<(N/2+1); i++)
      fprintf(outFile, "%E\n", psd[i]/(double)(nSegs));
 
   /* Housekeeping */
   
   printf("\nDone!\n"); 
   fclose(inFile);
   fclose(outFile);
 
   fftw_destroy_plan(plan);
   fftw_free(in); 
   fftw_free(out);
   free(psd);
   }

---

thank you very much thylacine1975.

that is what im looking for and im gonna try this library. it looks perfect, awesomely understandable for me
thanks for helping me! take care~

 

[mashhur]

---------- Post added at 02:53 ---------- Previous post was at 02:51 ----------

If you only need to know the signal's strength in DB, then you don't need an FFT or Hilbert transform.

Can't you just import the data into a spreadsheet, calculate the average value, subtract that from all the datapoints to get the AC, then calculate the RMS value and convert to dB?

hi godfreyl,

im curious to know your way. could you please tell me how to implement this step by step.
if i have 8192 values in 0.12us and sampling frequency is 1MHz, how do i get signal's RMS?

1. What is a Spreadsheet?
2. Average value (1+2+ .... + N)/NumberOfParameters (in my case NumberOfParameters=8192)
3. to get AC? what does AC stand for?
4. RMS = 10lg(AC)?
5. How to convert RMS into dB?

thank you very much!

 

[mgate]

mashhur,
1- A spreadsheet is an Microsoft excel worksheet

2 - if you have a wave and subtract to it the average value of the same wave, you will get a wave with zero average value... thats probably the AC (i think)

3 - RMS stands for root medium square. and is the square root of the medium of the sum of the square of each parameter

anyway, I don't think this what you are looking for, but if you want to use excel to calculate the DB, please feel free to check this workbook available at:

version excel 2003
**broken link removed**

version excel 2007
**broken link removed**

check excel addin configuration and sheet's instructions at
**broken link removed**