kopia lustrzana https://github.com/Dsplib/libdspl-2.0
added new fft function (doesn't work)
Dsplib
rodzic
6f01ae80cd
commit
8378a3517d
|
|
@ -0,0 +1,24 @@
|
|||
#ifndef DSPL_INTERNAL_H
|
||||
#define DSPL_INTERNAL_H
|
||||
|
||||
|
||||
#define DSPL_FARROW_LAGRANGE_COEFF 0.16666666666666666666666666666667
|
||||
|
||||
#define DSPL_RAND_MOD_X1 2147483647
|
||||
#define DSPL_RAND_MOD_X2 2145483479
|
||||
|
||||
|
||||
void transpose(double* a, int n, int m, double* b);
|
||||
void transpose_cmplx(complex_t* a, int n, int m, complex_t* b);
|
||||
void transpose_hermite(complex_t* a, int n, int m, complex_t* b);
|
||||
|
||||
int fft_bit_reverse(complex_t* x, complex_t* y, int n, int p2);
|
||||
int fft_dit(fft_t *pfft, int n, complex_t* y);
|
||||
void fft_dit_krn(complex_t *x0, complex_t *x1, complex_t *w, int n,
|
||||
complex_t *y0, complex_t *y1);
|
||||
|
||||
int fft_p2(int n);
|
||||
void dft2 (complex_t *x, complex_t* y);
|
||||
|
||||
#endif
|
||||
|
||||
123
dspl/src/fft.c
123
dspl/src/fft.c
|
|
@ -21,18 +21,8 @@
|
|||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include "dspl.h"
|
||||
#include "dspl_internal.h"
|
||||
|
||||
int fft_bit_reverse(complex_t* x, complex_t* y, int n, int p2);
|
||||
|
||||
int fft_dit(fft_t *pfft, int n, complex_t* y);
|
||||
|
||||
void fft_dit_krn(complex_t *x0, complex_t *x1, complex_t *w, int n,
|
||||
complex_t *y0, complex_t *y1);
|
||||
|
||||
int fft_p2(int n);
|
||||
|
||||
|
||||
void dft2 (complex_t *x, complex_t* y)
|
||||
|
||||
|
||||
|
||||
|
|
@ -256,6 +246,116 @@ int DSPL_API fft_create(fft_t *pfft, int n)
|
|||
|
||||
|
||||
|
||||
|
||||
/*******************************************************************************
|
||||
composite FFT kernel
|
||||
*******************************************************************************/
|
||||
int DSPL_API fftn_krn(complex_t* t0, complex_t* t1, fft_t* p, int n, int addr)
|
||||
{
|
||||
int n1, n2, s, k;
|
||||
s = n;
|
||||
addr = 0;
|
||||
|
||||
n2 = 1;
|
||||
if(s%2 == 0)
|
||||
{
|
||||
n2 = 2;
|
||||
goto label_size;
|
||||
}
|
||||
|
||||
label_size:
|
||||
if(n2 == 1)
|
||||
return ERROR_FFT_SIZE;
|
||||
n1 = s / n2;
|
||||
memcpy(t1, t0, n*sizeof(complex_t));
|
||||
|
||||
transpose_cmplx(t1, n1, n2, t0);
|
||||
|
||||
if(n2 == 2)
|
||||
{
|
||||
for(k = 0; k < n1; k++)
|
||||
{
|
||||
dft2(t0+2*k, t1+2*k);
|
||||
}
|
||||
}
|
||||
|
||||
if(n1 > 1)
|
||||
{
|
||||
for(k =0; k < n; k++)
|
||||
{
|
||||
RE(t0[k]) = CMRE(t1[k], p->w[addr+k]);
|
||||
IM(t0[k]) = CMIM(t1[k], p->w[addr+k]);
|
||||
}
|
||||
transpose_cmplx(t0, n2, n1, t1);
|
||||
|
||||
for(k = 0; k < n2; k++)
|
||||
{
|
||||
fftn_krn(t1, t0, p, s, addr+n);
|
||||
}
|
||||
transpose_cmplx(t0, n1, n2, t1);
|
||||
}
|
||||
return RES_OK;
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/*******************************************************************************
|
||||
FFT create for composite N
|
||||
*******************************************************************************/
|
||||
int DSPL_API fftn_create(fft_t *pfft, int n)
|
||||
{
|
||||
|
||||
int n1, n2, addr, s, k, m, nw;
|
||||
double phi;
|
||||
s = n;
|
||||
nw = addr = 0;
|
||||
|
||||
while(s > 1)
|
||||
{
|
||||
n2 = 1;
|
||||
if(s%2 == 0)
|
||||
{
|
||||
n2 = 2;
|
||||
goto label_size;
|
||||
}
|
||||
|
||||
|
||||
label_size:
|
||||
if(n2 == 1)
|
||||
return ERROR_FFT_SIZE;
|
||||
n1 = s / n2;
|
||||
nw += s;
|
||||
pfft->w = pfft->w ? (complex_t*) realloc(pfft->w, nw*sizeof(complex_t)):
|
||||
(complex_t*) malloc( nw*sizeof(complex_t));
|
||||
|
||||
for(k = 0; k < n2; k++)
|
||||
{
|
||||
for(m = 0; m < n1; m++)
|
||||
{
|
||||
phi = - M_2PI * (double)(k*m) / (double)s;
|
||||
RE(pfft->w[addr]) = cos(phi);
|
||||
IM(pfft->w[addr]) = sin(phi);
|
||||
addr++;
|
||||
}
|
||||
}
|
||||
s /= n2;
|
||||
}
|
||||
|
||||
pfft->t0 = pfft->t0 ? (complex_t*) realloc(pfft->t0, n*sizeof(complex_t)):
|
||||
(complex_t*) malloc( n*sizeof(complex_t));
|
||||
|
||||
|
||||
pfft->t1 = pfft->t1 ? (complex_t*) realloc(pfft->t1, n*sizeof(complex_t)):
|
||||
(complex_t*) malloc( n*sizeof(complex_t));
|
||||
|
||||
return RES_OK;
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
/*******************************************************************************
|
||||
FFT decimation in time
|
||||
*******************************************************************************/
|
||||
|
|
@ -347,6 +447,7 @@ void DSPL_API fft_free(fft_t *pfft)
|
|||
free(pfft->t0);
|
||||
if(pfft->t1)
|
||||
free(pfft->t1);
|
||||
memset(pfft, 0, sizeof(fft_t));
|
||||
}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -21,12 +21,9 @@
|
|||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include "dspl.h"
|
||||
#include "dspl_internal.h"
|
||||
|
||||
|
||||
void transpose(double* a, int n, int m, double* b);
|
||||
void transpose_cmplx(complex_t* a, int n, int m, complex_t* b);
|
||||
void transpose_hermite(complex_t* a, int n, int m, complex_t* b);
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -24,10 +24,9 @@
|
|||
#include <time.h>
|
||||
|
||||
#include "dspl.h"
|
||||
#include "dspl_internal.h"
|
||||
|
||||
|
||||
#define DSPL_RAND_MOD_X1 2147483647
|
||||
#define DSPL_RAND_MOD_X2 2145483479
|
||||
|
||||
|
||||
/******************************************************************************
|
||||
|
|
|
|||
|
|
@ -24,8 +24,8 @@
|
|||
#include <string.h>
|
||||
#include <math.h>
|
||||
#include "dspl.h"
|
||||
#include "dspl_internal.h"
|
||||
|
||||
#define DSPL_FARROW_LAGRANGE_COEFF 0.16666666666666666666666666666667
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -42,7 +42,7 @@ p_butter_ap_zp butter_ap_zp ;
|
|||
p_cheby_poly1 cheby_poly1 ;
|
||||
p_cheby_poly2 cheby_poly2 ;
|
||||
p_cheby1_ap cheby1_ap ;
|
||||
p_cheby1_ap_zp cheby1_ap_zp ;
|
||||
p_cheby1_ap_zp cheby1_ap_zp ;
|
||||
p_cheby2_ap cheby2_ap ;
|
||||
p_cheby2_ap_zp cheby2_ap_zp ;
|
||||
p_cmplx2re cmplx2re ;
|
||||
|
|
@ -82,6 +82,8 @@ p_fft_create fft_create ;
|
|||
p_fft_free fft_free ;
|
||||
p_fft_shift fft_shift ;
|
||||
p_fft_shift_cmplx fft_shift_cmplx ;
|
||||
p_fftn_create fftn_create ;
|
||||
p_fftn_krn fftn_krn ;
|
||||
p_flipip flipip ;
|
||||
p_flipip_cmplx flipip_cmplx ;
|
||||
p_fourier_series_dec fourier_series_dec ;
|
||||
|
|
@ -126,6 +128,8 @@ p_writetxt_3d writetxt_3d ;
|
|||
p_writetxt_3dline writetxt_3dline ;
|
||||
p_writetxt_cmplx_im writetxt_cmplx_im ;
|
||||
p_writetxt_cmplx_re writetxt_cmplx_re ;
|
||||
|
||||
|
||||
#endif //BUILD_LIB
|
||||
|
||||
|
||||
|
|
@ -223,6 +227,8 @@ void* dspl_load()
|
|||
LOAD_FUNC(fft_free);
|
||||
LOAD_FUNC(fft_shift);
|
||||
LOAD_FUNC(fft_shift_cmplx);
|
||||
LOAD_FUNC(fftn_create);
|
||||
LOAD_FUNC(fftn_krn);
|
||||
LOAD_FUNC(flipip);
|
||||
LOAD_FUNC(flipip_cmplx);
|
||||
LOAD_FUNC(fourier_series_dec);
|
||||
|
|
|
|||
|
|
@ -317,10 +317,10 @@ DECLARE_FUNC(int, dft_cmplx, complex_t*
|
|||
COMMA int
|
||||
COMMA complex_t*);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(double, dmod, double
|
||||
DECLARE_FUNC(double, dmod, double
|
||||
COMMA double);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(void, dspl_info, void);
|
||||
DECLARE_FUNC(void, dspl_info, void);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(int, ellip_acd, double* w
|
||||
COMMA int n
|
||||
|
|
@ -444,7 +444,7 @@ DECLARE_FUNC(int, fft_cmplx, complex_t*
|
|||
DECLARE_FUNC(int, fft_create, fft_t*
|
||||
COMMA int);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(void, fft_free, fft_t*);
|
||||
DECLARE_FUNC(void, fft_free, fft_t*);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(int, fft_shift, double*
|
||||
COMMA int n
|
||||
|
|
@ -454,6 +454,15 @@ DECLARE_FUNC(int, fft_shift_cmplx, complex_t*
|
|||
COMMA int
|
||||
COMMA complex_t*);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(int, fftn_create, fft_t* pfft
|
||||
COMMA int n);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(int, fftn_krn, complex_t* t0
|
||||
COMMA complex_t* t1
|
||||
COMMA fft_t* p
|
||||
COMMA int n
|
||||
COMMA int addr);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(int, flipip, double*
|
||||
COMMA int);
|
||||
//------------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -42,7 +42,7 @@ p_butter_ap_zp butter_ap_zp ;
|
|||
p_cheby_poly1 cheby_poly1 ;
|
||||
p_cheby_poly2 cheby_poly2 ;
|
||||
p_cheby1_ap cheby1_ap ;
|
||||
p_cheby1_ap_zp cheby1_ap_zp ;
|
||||
p_cheby1_ap_zp cheby1_ap_zp ;
|
||||
p_cheby2_ap cheby2_ap ;
|
||||
p_cheby2_ap_zp cheby2_ap_zp ;
|
||||
p_cmplx2re cmplx2re ;
|
||||
|
|
@ -82,6 +82,8 @@ p_fft_create fft_create ;
|
|||
p_fft_free fft_free ;
|
||||
p_fft_shift fft_shift ;
|
||||
p_fft_shift_cmplx fft_shift_cmplx ;
|
||||
p_fftn_create fftn_create ;
|
||||
p_fftn_krn fftn_krn ;
|
||||
p_flipip flipip ;
|
||||
p_flipip_cmplx flipip_cmplx ;
|
||||
p_fourier_series_dec fourier_series_dec ;
|
||||
|
|
@ -126,6 +128,8 @@ p_writetxt_3d writetxt_3d ;
|
|||
p_writetxt_3dline writetxt_3dline ;
|
||||
p_writetxt_cmplx_im writetxt_cmplx_im ;
|
||||
p_writetxt_cmplx_re writetxt_cmplx_re ;
|
||||
|
||||
|
||||
#endif //BUILD_LIB
|
||||
|
||||
|
||||
|
|
@ -223,6 +227,8 @@ void* dspl_load()
|
|||
LOAD_FUNC(fft_free);
|
||||
LOAD_FUNC(fft_shift);
|
||||
LOAD_FUNC(fft_shift_cmplx);
|
||||
LOAD_FUNC(fftn_create);
|
||||
LOAD_FUNC(fftn_krn);
|
||||
LOAD_FUNC(flipip);
|
||||
LOAD_FUNC(flipip_cmplx);
|
||||
LOAD_FUNC(fourier_series_dec);
|
||||
|
|
|
|||
|
|
@ -317,10 +317,10 @@ DECLARE_FUNC(int, dft_cmplx, complex_t*
|
|||
COMMA int
|
||||
COMMA complex_t*);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(double, dmod, double
|
||||
DECLARE_FUNC(double, dmod, double
|
||||
COMMA double);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(void, dspl_info, void);
|
||||
DECLARE_FUNC(void, dspl_info, void);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(int, ellip_acd, double* w
|
||||
COMMA int n
|
||||
|
|
@ -444,7 +444,7 @@ DECLARE_FUNC(int, fft_cmplx, complex_t*
|
|||
DECLARE_FUNC(int, fft_create, fft_t*
|
||||
COMMA int);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(void, fft_free, fft_t*);
|
||||
DECLARE_FUNC(void, fft_free, fft_t*);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(int, fft_shift, double*
|
||||
COMMA int n
|
||||
|
|
@ -454,6 +454,15 @@ DECLARE_FUNC(int, fft_shift_cmplx, complex_t*
|
|||
COMMA int
|
||||
COMMA complex_t*);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(int, fftn_create, fft_t* pfft
|
||||
COMMA int n);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(int, fftn_krn, complex_t* t0
|
||||
COMMA complex_t* t1
|
||||
COMMA fft_t* p
|
||||
COMMA int n
|
||||
COMMA int addr);
|
||||
//------------------------------------------------------------------------------
|
||||
DECLARE_FUNC(int, flipip, double*
|
||||
COMMA int);
|
||||
//------------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -0,0 +1,33 @@
|
|||
clear all;
|
||||
close all;
|
||||
clc;
|
||||
|
||||
N = 6;
|
||||
N1 = 3;
|
||||
N2 = 2;
|
||||
|
||||
% входной сигнал это вектор столбец размерности [N x 1]
|
||||
x = (0:N-1)';
|
||||
|
||||
for n1 = 0 : N1-1
|
||||
for k2 = 0 : N2-1
|
||||
W(n1 + 1, k2 + 1) = exp(-2i * pi * n1 * k2 / N);
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
A = reshape(x, N2, N1);
|
||||
|
||||
B = A.';
|
||||
D = fft(B);
|
||||
F = D.*W;
|
||||
G = F.';
|
||||
H = fft(G);
|
||||
P = H.';
|
||||
|
||||
y = [reshape(P,N,1), fft(x)]
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -0,0 +1,37 @@
|
|||
#include <stdio.h>
|
||||
#include <stdlib.h>
|
||||
#include <string.h>
|
||||
#include "dspl.h"
|
||||
|
||||
#define N 16
|
||||
int main()
|
||||
{
|
||||
void* handle; // DSPL handle
|
||||
handle = dspl_load(); // Load DSPL function
|
||||
int n;
|
||||
|
||||
complex_t x[N];
|
||||
complex_t y[N]; // DFT
|
||||
fft_t pfft; // FFT object
|
||||
|
||||
//clear fft object
|
||||
memset(&pfft, 0, sizeof(fft_t));
|
||||
|
||||
// Create FFT object
|
||||
fftn_create(&pfft, N);
|
||||
|
||||
for(n = 0; n < N; n++)
|
||||
{
|
||||
RE(x[n]) = (double)n;
|
||||
IM(x[n]) = 0.0;
|
||||
}
|
||||
|
||||
fftn_krn(x, y, &pfft, N, 0);
|
||||
|
||||
|
||||
fft_free(&pfft); // clear FFT object
|
||||
|
||||
dspl_free(handle); // free dspl handle
|
||||
return 0;
|
||||
}
|
||||
|
||||
Ładowanie…
Reference in New Issue