RTXI  2.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
welcpdgm.cpp
Go to the documentation of this file.
1 //
2 // File = welcpdgm.cpp
3 //
4 
5 #include "welcpdgm.h"
6 #include "complex.h"
7 #include "fft.h"
8 #include "misdefs.h"
9 #include <stdlib.h>
10 #include <string.h>
11 
13  double samp_intvl, int num_samps_per_seg,
14  int shift_between_segs, int fft_len,
15  GenericWindow* data_wind,
16  int num_segs_to_avg)
17  : PsdEstimate(num_samps_per_seg, samp_intvl)
18 {
19  int samp_idx, seg_idx, ovrlap_len;
20  complex *time_seg, *ovrlap_seg, *new_seg, *freq_seg;
21  complex* windowed_seg;
22  double* win_seq;
23  double scale_factor;
24 
25  std::cout << "in WelchPeriodogram" << std::endl;
26 
27  Psd_Est = new double[fft_len];
28  time_seg = new complex[num_samps_per_seg];
29  windowed_seg = new complex[num_samps_per_seg];
30  freq_seg = new complex[fft_len];
31 
32  if ((time_seg == NULL) || (windowed_seg == NULL) || (freq_seg == NULL)) {
33  std::cout << "Allocation error in freq_seg in WelchPeriodogram"
34  << std::endl;
35  exit(99);
36  }
37 
38  for (samp_idx = 0; samp_idx < fft_len; samp_idx++) {
39  Psd_Est[samp_idx] = 0.0;
40  }
41  if (data_wind != NULL)
42  win_seq = data_wind->GetDataWindow();
43 
44  ovrlap_seg = time_seg + shift_between_segs;
45  ovrlap_len = num_samps_per_seg - shift_between_segs;
46  new_seg = time_seg + ovrlap_len;
47 
48  // prefill overlap segment of buffer
49  signal_source->GetNextSegment(ovrlap_seg, ovrlap_len);
50  for (seg_idx = 0; seg_idx < num_segs_to_avg; seg_idx++) {
51  // copy overlap samples down to start of buffer
52  memmove(time_seg, ovrlap_seg,
53  (sizeof(complex) / sizeof(char)) * ovrlap_len);
54 
55  // get new samples
56  signal_source->GetNextSegment(new_seg, shift_between_segs);
57 
58  // apply window sequence to segment of time signal
59  if (data_wind == NULL) {
60  for (samp_idx = 0; samp_idx < num_samps_per_seg; samp_idx++) {
61  windowed_seg[samp_idx] = time_seg[samp_idx];
62  }
63  } else {
64  for (samp_idx = 0; samp_idx < num_samps_per_seg; samp_idx++) {
65  windowed_seg[samp_idx] = win_seq[samp_idx] * time_seg[samp_idx];
66  }
67  }
68 
69  //-----------------------------------
70  // compute FFT of windowed segment
71  fft(windowed_seg, freq_seg, num_samps_per_seg, fft_len);
72 
73  for (samp_idx = 0; samp_idx < fft_len; samp_idx++) {
74  Psd_Est[samp_idx] += mag_sqrd(freq_seg[samp_idx]);
75  }
76  }
77  scale_factor = (float)num_segs_to_avg * num_samps_per_seg / samp_intvl;
78  for (samp_idx = 0; samp_idx < fft_len; samp_idx++) {
79  Psd_Est[samp_idx] /= scale_factor;
80  }
81  delete[] time_seg;
82  delete[] freq_seg;
83  delete[] windowed_seg;
84  return;
85 }
86 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
double * GetDataWindow(void)
Definition: gen_win.cpp:74
WelchPeriodogram(SignalSource *signal_source, double samp_intvl, int num_samps_per_seg, int shift_between_segs, int fft_len, GenericWindow *data_wind, int num_segs_to_avg)
Definition: welcpdgm.cpp:12
void fft(complex *time_signal, complex *sample_spectrum, int num_samps)
Definition: fft.cpp:24
double mag_sqrd(const complex _z)
Definition: complex.h:128
double * Psd_Est
Definition: psd_est.h:25
virtual void GetNextSegment(complex *, int)
Definition: sig_src.h:19