ASCEE/lasp
ASCEE/lasp
5
1
Fork 0
Code Issues 5 Pull Requests Actions Releases 14 Activity

Fixed scaling issue in FFT

Browse Source
This commit is contained in:
Casper Jansen 2021-08-13 11:40:25 +02:00
parent 4c0f399505
commit c98b10b83a
1 changed files with 16 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
lasp/c/lasp_ps.c
View File

@ -1,7 +1,7 @@
// lasp_ps.c // lasp_ps.c
// //
// Author: J.A. de Jong -ASCEE // Author: J.A. de Jong -ASCEE
// //
// Description: // Description:
// //
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
@ -21,7 +21,7 @@ typedef struct PowerSpectra_s {
PowerSpectra* PowerSpectra_alloc(const us nfft, PowerSpectra* PowerSpectra_alloc(const us nfft,
const WindowType wt) { const WindowType wt) {
fsTRACE(15); fsTRACE(15);
int rv; int rv;
@ -48,7 +48,7 @@ PowerSpectra* PowerSpectra_alloc(const us nfft,
/* Allocate vectors and matrices */ /* Allocate vectors and matrices */
ps->window = vd_alloc(nfft); ps->window = vd_alloc(nfft);
rv = window_create(wt,&ps->window,&ps->win_pow); rv = window_create(wt,&ps->window,&ps->win_pow);
check_overflow_vx(ps->window); check_overflow_vx(ps->window);
if(rv!=0) { if(rv!=0) {
@ -70,11 +70,11 @@ void PowerSpectra_free(PowerSpectra* ps) {
void PowerSpectra_compute(const PowerSpectra* ps, void PowerSpectra_compute(const PowerSpectra* ps,
const dmat * timedata, const dmat * timedata,
cmat * result) { cmat * result) {
fsTRACE(15); fsTRACE(15);
dbgassert(ps && timedata && result,NULLPTRDEREF); dbgassert(ps && timedata && result,NULLPTRDEREF);
const us nchannels = timedata->n_cols; const us nchannels = timedata->n_cols;
const us nfft = Fft_nfft(ps->fft); const us nfft = Fft_nfft(ps->fft);
uVARTRACE(15,nchannels); uVARTRACE(15,nchannels);
@ -100,13 +100,13 @@ void PowerSpectra_compute(const PowerSpectra* ps,
vd column_work = dmat_column(&timedata_work,i); vd column_work = dmat_column(&timedata_work,i);
vd_elem_prod(&column_work,&column,&ps->window); vd_elem_prod(&column_work,&column,&ps->window);
vd_free(&column); vd_free(&column);
vd_free(&column_work); vd_free(&column_work);
} }
check_overflow_xmat(timedata_work); check_overflow_xmat(timedata_work);
/* print_dmat(&timedata_work); */ /* print_dmat(&timedata_work); */
/* Compute fft of the time data */ /* Compute fft of the time data */
cmat fft_work = cmat_alloc(nfft/2+1,nchannels); cmat fft_work = cmat_alloc(nfft/2+1,nchannels);
Fft_fft(ps->fft, Fft_fft(ps->fft,
@ -114,19 +114,22 @@ void PowerSpectra_compute(const PowerSpectra* ps,
&fft_work); &fft_work);
dmat_free(&timedata_work); dmat_free(&timedata_work);
TRACE(15,"fft done"); TRACE(15,"fft done");
/* Scale fft such that power is easily computed */ /* Scale fft such that power is easily computed
const c scale_fac = d_sqrt(2/win_pow)/nfft; * - Multiply power spectral density by 2 except at f=0 and f=fNq
* - Divide by energy of window function = nfft * window_power
* - .. sqrt(factors) because it is applied to output fft instead of psd */
const c scale_fac = d_sqrt(2/(nfft*win_pow));
cmat_scale(&fft_work,scale_fac); cmat_scale(&fft_work,scale_fac);
TRACE(15,"scale done"); TRACE(15,"scale done");
for(us i=0;i< nchannels;i++) { for(us i=0;i< nchannels;i++) {
/* Multiply DC term with 1/sqrt(2) */ /* Multiply DC term by 1/sqrt(2) */
*getcmatval(&fft_work,0,i) *= 1/d_sqrt(2.)+0*I; *getcmatval(&fft_work,0,i) *= 1/d_sqrt(2.)+0*I;
/* Multiply Nyquist term with 1/sqrt(2) */ /* Multiply Nyquist term by 1/sqrt(2) */
*getcmatval(&fft_work,nfft/2,i) *= 1/d_sqrt(2.)+0*I; *getcmatval(&fft_work,nfft/2,i) *= 1/d_sqrt(2.)+0*I;
} }
check_overflow_xmat(fft_work); check_overflow_xmat(fft_work);