File indexing completed on 2024-05-12 15:27:06
0001 /*************************************************************************** 0002 File : nsl_dft.c 0003 Project : LabPlot 0004 Description : NSL discrete Fourier transform functions 0005 -------------------------------------------------------------------- 0006 Copyright : (C) 2016 by Stefan Gerlach (stefan.gerlach@uni.kn) 0007 0008 ***************************************************************************/ 0009 0010 /*************************************************************************** 0011 * * 0012 * This program is free software; you can redistribute it and/or modify * 0013 * it under the terms of the GNU General Public License as published by * 0014 * the Free Software Foundation; either version 2 of the License, or * 0015 * (at your option) any later version. * 0016 * * 0017 * This program is distributed in the hope that it will be useful, * 0018 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 0019 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 0020 * GNU General Public License for more details. * 0021 * * 0022 * You should have received a copy of the GNU General Public License * 0023 * along with this program; if not, write to the Free Software * 0024 * Foundation, Inc., 51 Franklin Street, Fifth Floor, * 0025 * Boston, MA 02110-1301 USA * 0026 * * 0027 ***************************************************************************/ 0028 0029 #include "nsl_dft.h" 0030 #include "nsl_common.h" 0031 #include <gsl/gsl_fft_real.h> 0032 #include <gsl/gsl_fft_halfcomplex.h> 0033 #ifdef HAVE_FFTW3 0034 #include <fftw3.h> 0035 #endif 0036 0037 const char* nsl_dft_result_type_name[] = {i18n("Magnitude"), i18n("Amplitude"), i18n("real part"), i18n("imaginary part"), i18n("Power"), i18n("Phase"), 0038 i18n("Amplitude in dB"), i18n("normalized amplitude in dB"), i18n("Magnitude squared"), i18n("Amplitude squared"), i18n("raw")}; 0039 const char* nsl_dft_xscale_name[] = {i18n("Frequency"), i18n("Index"), i18n("Period")}; 0040 0041 int nsl_dft_transform_window(double data[], size_t stride, size_t n, int two_sided, nsl_dft_result_type type, nsl_sf_window_type window_type) { 0042 /* apply window function */ 0043 if (window_type != nsl_sf_window_uniform) 0044 nsl_sf_apply_window(data, n, window_type); 0045 0046 /* transform */ 0047 int status = nsl_dft_transform(data, stride, n, two_sided, type); 0048 0049 return status; 0050 } 0051 0052 int nsl_dft_transform(double data[], size_t stride, size_t n, int two_sided, nsl_dft_result_type type) { 0053 if (n < 2) // we need at least 2 points 0054 return 1; 0055 size_t i; 0056 double* result = (double*)malloc(2*n*sizeof(double)); 0057 size_t N = n/2; /* number of resulting data points */ 0058 if(two_sided) 0059 N = n; 0060 #ifdef HAVE_FFTW3 0061 /* stride ignored */ 0062 (void)stride; 0063 0064 fftw_plan plan = fftw_plan_dft_r2c_1d(n, data, (fftw_complex *) result, FFTW_ESTIMATE); 0065 fftw_execute(plan); 0066 fftw_destroy_plan(plan); 0067 0068 /* 2. unpack data */ 0069 if(two_sided) { 0070 for(i = 1; i < n-i; i++) { 0071 result[2*(n - i)] = result[2*i]; 0072 result[2*(n - i)+1] = -result[2*i+1]; 0073 } 0074 if (i == n - i) { 0075 result[2*i] = result[2*(n-1)]; 0076 result[2*i+1] = 0; 0077 } 0078 } 0079 #else 0080 /* 1. transform */ 0081 gsl_fft_real_wavetable *real = gsl_fft_real_wavetable_alloc(n); 0082 gsl_fft_real_workspace *work = gsl_fft_real_workspace_alloc(n); 0083 0084 gsl_fft_real_transform(data, stride, n, real, work); 0085 gsl_fft_real_wavetable_free(real); 0086 gsl_fft_real_workspace_free(work); 0087 0088 /* 2. unpack data */ 0089 gsl_fft_halfcomplex_unpack(data, result, stride, n); 0090 #endif 0091 0092 /* 3. write result */ 0093 switch(type) { 0094 case nsl_dft_result_magnitude: 0095 for (i = 0; i < N; i++) 0096 data[i] = sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1])); 0097 break; 0098 case nsl_dft_result_amplitude: 0099 for (i = 0; i < N; i++) { 0100 data[i] = sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/n; 0101 if(i>0) 0102 data[i] *= 2.; 0103 } 0104 break; 0105 case nsl_dft_result_real: 0106 for (i = 0; i < N; i++) 0107 data[i] = result[2*i]; 0108 break; 0109 case nsl_dft_result_imag: 0110 for (i = 0; i < N; i++) 0111 data[i] = result[2*i+1]; 0112 break; 0113 case nsl_dft_result_power: 0114 for (i = 0; i < N; i++) { 0115 data[i] = (gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/n; 0116 if (i > 0) 0117 data[i] *= 2.; 0118 } 0119 break; 0120 case nsl_dft_result_phase: 0121 for (i = 0; i < N; i++) 0122 data[i] = -atan2(result[2*i+1],result[2*i]); 0123 break; 0124 case nsl_dft_result_dB: 0125 for (i = 0; i < N; i++) { 0126 if (i == 0) 0127 data[i] = 20.*log10(sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/(double)n); 0128 else 0129 data[i] = 20.*log10(2.*sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/(double)n); 0130 } 0131 break; 0132 case nsl_dft_result_normdB: { 0133 double maxdB=0; 0134 for (i = 0; i < N; i++) { 0135 if (i == 0) 0136 data[i] = 20.*log10(sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/(double)n); 0137 else 0138 data[i] = 20.*log10(2.*sqrt(gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]))/(double)n); 0139 0140 if (i == 0 || maxdB < data[i]) 0141 maxdB = data[i]; 0142 } 0143 for (i = 0; i < N; i++) 0144 data[i] -= maxdB; 0145 break; 0146 } 0147 case nsl_dft_result_squaremagnitude: 0148 for (i = 0; i < N; i++) 0149 data[i] = gsl_pow_2(result[2*i])+gsl_pow_2(result[2*i+1]); 0150 break; 0151 case nsl_dft_result_squareamplitude: 0152 for (i = 0; i < N; i++) { 0153 data[i] = (gsl_pow_2(result[2*i]) + gsl_pow_2(result[2*i+1]))/gsl_pow_2((double)n); 0154 if (i > 0) 0155 data[i] *= 4.; 0156 } 0157 break; 0158 case nsl_dft_result_raw: 0159 break; 0160 } 0161 0162 free(result); 0163 0164 return 0; 0165 } 0166