File indexing completed on 2024-05-12 03:47:52
0001 /* 0002 File : nsl_hilbert.c 0003 Project : LabPlot 0004 Description : NSL Hilbert transform 0005 -------------------------------------------------------------------- 0006 SPDX-FileCopyrightText: 2021 Stefan Gerlach <stefan.gerlach@uni.kn> 0007 SPDX-License-Identifier: GPL-2.0-or-later 0008 */ 0009 0010 #include "nsl_hilbert.h" 0011 #include "nsl_common.h" 0012 #include <gsl/gsl_fft_complex.h> 0013 #include <gsl/gsl_fft_halfcomplex.h> 0014 #ifdef HAVE_FFTW3 0015 #include <fftw3.h> 0016 #endif 0017 0018 const char* nsl_hilbert_result_type_name[] = {i18n("Imaginary Part"), i18n("Envelope")}; 0019 0020 /* algorithm from https://de.wikipedia.org/wiki/Hilbert-Transformation#Berechnung_%C3%BCber_Fouriertransformation */ 0021 int nsl_hilbert_transform(double data[], size_t stride, size_t n, nsl_hilbert_result_type type) { 0022 if (n < 2) /* we need at least 2 points */ 0023 return 1; 0024 0025 /* 1. DFT of data: dft_transform returns gsl_halfcomplex (raw) */ 0026 nsl_dft_transform(data, stride, n, 1, nsl_dft_result_raw); 0027 0028 const size_t N = 2 * n; 0029 double* result = (double*)malloc(N * sizeof(double)); 0030 gsl_fft_halfcomplex_unpack(data, result, stride, n); 0031 0032 size_t i; 0033 /* for (i = 0; i < N; i++) 0034 printf(" %g", result[i]); 0035 printf("\n"); 0036 */ 0037 /* 2. manipulate values */ 0038 for (i = 2; i < 2 * ceil(n / 2.); i++) 0039 result[i] *= 2; 0040 0041 /* for (i = 0; i < N; i++) 0042 printf(" %g", result[i]); 0043 printf("\n"); 0044 */ 0045 for (i = n + 1; i < N; i++) 0046 result[i] = 0; 0047 0048 /* for (i = 0; i < N; i++) 0049 printf(" %g", result[i]); 0050 printf("\n"); 0051 */ 0052 /* 3. back transform */ 0053 #ifdef HAVE_FFTW3 0054 fftw_complex* o = (fftw_complex*)malloc(N * sizeof(double)); 0055 fftw_plan pb = fftw_plan_dft_1d(n, o, o, FFTW_BACKWARD, FFTW_ESTIMATE); 0056 0057 fftw_execute_dft(pb, (fftw_complex*)result, (fftw_complex*)result); 0058 fftw_destroy_plan(pb); 0059 free(o); 0060 #else 0061 gsl_fft_complex_workspace* work = gsl_fft_complex_workspace_alloc(n); 0062 gsl_fft_complex_wavetable* hc = gsl_fft_complex_wavetable_alloc(n); 0063 gsl_fft_complex_inverse(result, 1, n, hc, work); 0064 gsl_fft_complex_wavetable_free(hc); 0065 gsl_fft_complex_workspace_free(work); 0066 #endif 0067 0068 /* for (i = 0; i < N; i++) 0069 #ifdef HAVE_FFTW3 0070 printf(" %g", result[i]/n); 0071 #else 0072 printf(" %g", result[i]); 0073 #endif 0074 printf("\n"); 0075 */ 0076 /* copy back (Hilbert transform is imag part) */ 0077 switch (type) { 0078 case nsl_hilbert_result_imag: 0079 for (i = 0; i < n; i++) 0080 #ifdef HAVE_FFTW3 0081 data[i] = result[2 * i + 1] / n; 0082 #else 0083 data[i] = result[2 * i + 1]; 0084 #endif 0085 break; 0086 case nsl_hilbert_result_envelope: 0087 for (i = 0; i < n; i++) 0088 #ifdef HAVE_FFTW3 0089 data[i] = hypot(result[2 * i], result[2 * i + 1]) / n; 0090 #else 0091 data[i] = hypot(result[2 * i], result[2 * i + 1]); 0092 #endif 0093 } 0094 0095 free(result); 0096 0097 return 0; 0098 }