mirror
/
libdspl-2.0
kopia lustrzana https://github.com/Dsplib/libdspl-2.0
1
0
Forkuj 0
Kod Zgłoszenia Projekty Wydania Wiki Aktywność

new fft and solved warnings

pull/6/merge
Dsplib 2018-12-11 23:58:50 +03:00
rodzic a80958c04b
commit 2a26f8aa4f
22 zmienionych plików z 392 dodań i 503 usunięć
Pokaż statystyki Ściągnij plik aktualizacji Ściągnij plik porównania Expand all files Collapse all files

14
dspl/src/array.c
Wyświetl plik

@ -162,6 +162,7 @@ int DSPL_API verif(double* x, double* y, size_t n, double eps, double* err)
{
double d, maxd;
size_t k;
int res;
if(!x || !y)
return ERROR_PTR;
if(n < 1)
@ -185,11 +186,11 @@ int DSPL_API verif(double* x, double* y, size_t n, double eps, double* err)
*err = maxd;
if(maxd > eps)
err = DSPL_VERIF_FAILED;
res = DSPL_VERIF_FAILED;
else
err = DSPL_VERIF_SUCCESS;
res = DSPL_VERIF_SUCCESS;
return err;
return res;
}
@ -204,6 +205,7 @@ int DSPL_API verif_cmplx(complex_t* x, complex_t* y, size_t n,
complex_t d;
double mx, md, maxd;
size_t k;
int res;
if(!x || !y)
return ERROR_PTR;
if(n < 1)
@ -230,9 +232,9 @@ int DSPL_API verif_cmplx(complex_t* x, complex_t* y, size_t n,
*err = maxd;
if(maxd > eps)
err = DSPL_VERIF_FAILED;
res = DSPL_VERIF_FAILED;
else
err = DSPL_VERIF_SUCCESS;
res = DSPL_VERIF_SUCCESS;
return err;
return res;
}

7
dspl/src/dspl_internal.h
Wyświetl plik

@ -12,12 +12,7 @@ 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);
int fft_krn(complex_t* t0, complex_t* t1, fft_t* p, int n, int addr);
void dft2 (complex_t *x, complex_t* y);

232
dspl/src/fft.c
Wyświetl plik

@ -46,7 +46,7 @@ int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
for(k = 0; k < n; k++)
IM(pfft->t1[k]) = -IM(pfft->t1[k]);
err = fft_dit(pfft, n, y);
err = fft_krn(pfft->t1, pfft->t0, pfft, n, 0);
if(err!=RES_OK)
return err;
@ -54,8 +54,8 @@ int DSPL_API ifft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
norm = 1.0 / (double)n;
for(k = 0; k < n; k++)
{
RE(y[k]) = RE(y[k])*norm;
IM(y[k]) = -IM(y[k])*norm;
RE(y[k]) = RE(pfft->t0[k])*norm;
IM(y[k]) = -IM(pfft->t0[k])*norm;
}
return RES_OK;
}
@ -80,7 +80,7 @@ int DSPL_API fft(double *x, int n, fft_t* pfft, complex_t* y)
re2cmplx(x, n, pfft->t1);
return fft_dit(pfft, n, y);
return fft_krn(pfft->t1, y, pfft, n, 0);
}
@ -105,135 +105,21 @@ int DSPL_API fft_cmplx(complex_t *x, int n, fft_t* pfft, complex_t* y)
memcpy(pfft->t1, x, n*sizeof(complex_t));
return fft_dit(pfft, n, y);
return fft_krn(pfft->t1, y, pfft, n, 0);
}
/*******************************************************************************
FFT bit reverse
*******************************************************************************/
int fft_bit_reverse(complex_t* x, complex_t* y, int n, int p2)
{
static unsigned char rb_table[256] =
{
0x00, 0x80, 0x40, 0xC0, 0x20, 0xA0, 0x60, 0xE0,
0x10, 0x90, 0x50, 0xD0, 0x30, 0xB0, 0x70, 0xF0,
0x08, 0x88, 0x48, 0xC8, 0x28, 0xA8, 0x68, 0xE8,
0x18, 0x98, 0x58, 0xD8, 0x38, 0xB8, 0x78, 0xF8,
0x04, 0x84, 0x44, 0xC4, 0x24, 0xA4, 0x64, 0xE4,
0x14, 0x94, 0x54, 0xD4, 0x34, 0xB4, 0x74, 0xF4,
0x0C, 0x8C, 0x4C, 0xCC, 0x2C, 0xAC, 0x6C, 0xEC,
0x1C, 0x9C, 0x5C, 0xDC, 0x3C, 0xBC, 0x7C, 0xFC,
0x02, 0x82, 0x42, 0xC2, 0x22, 0xA2, 0x62, 0xE2,
0x12, 0x92, 0x52, 0xD2, 0x32, 0xB2, 0x72, 0xF2,
0x0A, 0x8A, 0x4A, 0xCA, 0x2A, 0xAA, 0x6A, 0xEA,
0x1A, 0x9A, 0x5A, 0xDA, 0x3A, 0xBA, 0x7A, 0xFA,
0x06, 0x86, 0x46, 0xC6, 0x26, 0xA6, 0x66, 0xE6,
0x16, 0x96, 0x56, 0xD6, 0x36, 0xB6, 0x76, 0xF6,
0x0E, 0x8E, 0x4E, 0xCE, 0x2E, 0xAE, 0x6E, 0xEE,
0x1E, 0x9E, 0x5E, 0xDE, 0x3E, 0xBE, 0x7E, 0xFE,
0x01, 0x81, 0x41, 0xC1, 0x21, 0xA1, 0x61, 0xE1,
0x11, 0x91, 0x51, 0xD1, 0x31, 0xB1, 0x71, 0xF1,
0x09, 0x89, 0x49, 0xC9, 0x29, 0xA9, 0x69, 0xE9,
0x19, 0x99, 0x59, 0xD9, 0x39, 0xB9, 0x79, 0xF9,
0x05, 0x85, 0x45, 0xC5, 0x25, 0xA5, 0x65, 0xE5,
0x15, 0x95, 0x55, 0xD5, 0x35, 0xB5, 0x75, 0xF5,
0x0D, 0x8D, 0x4D, 0xCD, 0x2D, 0xAD, 0x6D, 0xED,
0x1D, 0x9D, 0x5D, 0xDD, 0x3D, 0xBD, 0x7D, 0xFD,
0x03, 0x83, 0x43, 0xC3, 0x23, 0xA3, 0x63, 0xE3,
0x13, 0x93, 0x53, 0xD3, 0x33, 0xB3, 0x73, 0xF3,
0x0B, 0x8B, 0x4B, 0xCB, 0x2B, 0xAB, 0x6B, 0xEB,
0x1B, 0x9B, 0x5B, 0xDB, 0x3B, 0xBB, 0x7B, 0xFB,
0x07, 0x87, 0x47, 0xC7, 0x27, 0xA7, 0x67, 0xE7,
0x17, 0x97, 0x57, 0xD7, 0x37, 0xB7, 0x77, 0xF7,
0x0F, 0x8F, 0x4F, 0xCF, 0x2F, 0xAF, 0x6F, 0xEF,
0x1F, 0x9F, 0x5F, 0xDF, 0x3F, 0xBF, 0x7F, 0xFF
};
if(!x || !y)
return ERROR_PTR;
if(n<1 || p2 < 1)
return ERROR_SIZE;
unsigned int v, c;
for(v = 0; v < (unsigned int)n; v++)
{
c = (unsigned int) ((rb_table[ v & 0xff] << 24) |
(rb_table[(v >> 8) & 0xff] << 16) |
(rb_table[(v >> 16) & 0xff] << 8) |
(rb_table[(v >> 24) & 0xff])) >> (32 - p2);
RE(y[c]) = RE(x[v]);
IM(y[c]) = IM(x[v]);
}
return RES_OK;
}
/*******************************************************************************
FFT create
*******************************************************************************/
int DSPL_API fft_create(fft_t *pfft, int n)
{
int p2, k,r,m,addr,s;
double phi;
p2 = fft_p2(n);
if(p2 < 1)
return ERROR_FFT_SIZE;
if(n < pfft->n+1)
return RES_OK;
pfft->n = n;
pfft->p2 = p2;
pfft->w = pfft->w ? (complex_t*) realloc(pfft->w, n*sizeof(complex_t)):
(complex_t*) malloc( n*sizeof(complex_t));
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));
m = 0;
addr = 0;
for(k = 0; k < p2; k++)
{
s = 1<<m;
for( r = 0; r < s; r++)
{
phi = -M_2PI *(double)r / (double)(2*s);
RE(pfft->w[addr+r]) = cos(phi);
IM(pfft->w[addr+r]) = sin(phi);
}
addr+=s;
m++;
}
return RES_OK;
}
/*******************************************************************************
composite FFT kernel
*******************************************************************************/
int DSPL_API fftn_krn(complex_t* t0, complex_t* t1, fft_t* p, int n, int addr)
int fft_krn(complex_t* t0, complex_t* t1, fft_t* p, int n, int addr)
{
int n1, n2, k;
complex_t *pw = p->w+addr;
n1 = 1;
if(n%16== 0) { n1 = 16; goto label_size; }
//if(n%16== 0) { n1 = 16; goto label_size; }
if(n%7 == 0) { n1 = 7; goto label_size; }
if(n%5 == 0) { n1 = 5; goto label_size; }
if(n%4 == 0) { n1 = 4; goto label_size; }
@ -289,7 +175,7 @@ label_size:
for(k = 0; k < n1; k++)
{
fftn_krn(t1+k*n2, t0+k*n2, p, n2, addr+n);
fft_krn(t1+k*n2, t0+k*n2, p, n2, addr+n);
}
transpose_cmplx(t0, n2, n1, t1);
}
@ -303,7 +189,7 @@ label_size:
/*******************************************************************************
FFT create for composite N
*******************************************************************************/
int DSPL_API fftn_create(fft_t *pfft, int n)
int DSPL_API fft_create(fft_t *pfft, int n)
{
int n1, n2, addr, s, k, m, nw;
@ -314,7 +200,7 @@ int DSPL_API fftn_create(fft_t *pfft, int n)
while(s > 1)
{
n2 = 1;
if(s%16== 0) { n2 = 16; goto label_size; }
//if(s%16== 0) { n2 = 16; goto label_size; }
if(s%7 == 0) { n2 = 7; goto label_size; }
if(s%5 == 0) { n2 = 5; goto label_size; }
if(s%4 == 0) { n2 = 4; goto label_size; }
@ -357,83 +243,6 @@ label_size:
/*******************************************************************************
FFT decimation in time
*******************************************************************************/
int fft_dit(fft_t *pfft, int n, complex_t* y)
{
int k,s,m,waddr, dm,p2;
complex_t *t, *t0, *t1, *w;
int err;
p2 = fft_p2(n);
if(p2<0)
return ERROR_FFT_SIZE;
t0 = pfft->t0;
t1 = pfft->t1;
w = pfft->w;
s = n>>1;
m = 1;
waddr = 0;
err = fft_bit_reverse(t1, t0, n, p2);
if(err!= RES_OK)
return err;
while(s)
{
dm = m<<1;
if(s > 1)
{
for(k = 0; k < n; k+=dm)
{
fft_dit_krn(t0+k, t0+k+m, w+waddr, m, t1+k, t1+k+m);
}
t = t1;
t1 = t0;
t0 = t;
waddr+=m;
m <<= 1;
}
else
{
fft_dit_krn(t0, t0+m, w+waddr, m, y, y+m);
}
s >>= 1;
}
return RES_OK;
}
/*******************************************************************************
FFT decimation in time kernel
*******************************************************************************/
void fft_dit_krn(complex_t *x0, complex_t *x1, complex_t *w, int n,
complex_t *y0, complex_t *y1)
{
int k;
complex_t mul;
for(k = 0; k < n; k++)
{
RE(mul) = CMRE(x1[k], w[k]);
IM(mul) = CMIM(x1[k], w[k]);
RE(y0[k]) = RE(x0[k]) + RE(mul);
IM(y0[k]) = IM(x0[k]) + IM(mul);
RE(y1[k]) = RE(x0[k]) - RE(mul);
IM(y1[k]) = IM(x0[k]) - IM(mul);
}
}
/*******************************************************************************
FFT free
@ -454,27 +263,6 @@ void DSPL_API fft_free(fft_t *pfft)
/*******************************************************************************
FFT power of 2
*******************************************************************************/
int fft_p2(int n)
{
int p = (n>0) ? n : 0;
int p2 = 0;
while((p = p>>1))
p2++;
if((1<<p2)!=n)
return -1;
return p2;
}
/*******************************************************************************
FFT shifting
*******************************************************************************/

4
dspl/src/fft_subkernel.c
Wyświetl plik

@ -218,8 +218,8 @@ void dft7 (complex_t *x, complex_t* y)
RE(sum[13]) = RE(sum[2]) + RE(sum[4]);
IM(sum[13]) = IM(sum[2]) + IM(sum[4]);
RE(sum[14]) = RE(sum[13]) + RE(x[6]);
IM(sum[14]) = IM(sum[13]) + IM(x[6]);
RE(sum[14]) = RE(sum[13]) + RE(sum[6]);
IM(sum[14]) = IM(sum[13]) + IM(sum[6]);
RE(sum[15]) = RE(sum[2]) - RE(sum[4]);
IM(sum[15]) = IM(sum[2]) - IM(sum[4]);

4
include/dspl.c
Wyświetl plik

@ -82,8 +82,6 @@ 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_integral_cmplx fourier_integral_cmplx ;
@ -231,8 +229,6 @@ 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_integral_cmplx);

9
include/dspl.h
Wyświetl plik

@ -457,15 +457,6 @@ 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);
//------------------------------------------------------------------------------

4
release/include/dspl.c
Wyświetl plik

@ -82,8 +82,6 @@ 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_integral_cmplx fourier_integral_cmplx ;
@ -231,8 +229,6 @@ 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_integral_cmplx);

9
release/include/dspl.h
Wyświetl plik

@ -457,15 +457,6 @@ 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);
//------------------------------------------------------------------------------

40
test/bin/gnuplot/fourier_transform_ex_gauss.plt 100755
Wyświetl plik

@ -0,0 +1,40 @@
set terminal wxt 0 size 560,480 enhanced font 'Verdana,8' position 0,0
unset key
set grid
set lmargin 8
set multiplot layout 2,1 rowsfirst
#
set xlabel 't, c'
set ylabel 's(t)'
plot[-10:10] 'dat/fourier_transform_ex_gauss_time_0.5.txt' with lines,\
'dat/fourier_transform_ex_gauss_time_1.0.txt' with lines,\
'dat/fourier_transform_ex_gauss_time_2.0.txt' with lines
#
set xlabel 'w, рад/c'
set ylabel 'S(w)'
plot[-4*pi:4*pi] 'dat/fourier_transform_ex_gauss_freq_0.5.txt' with lines,\
'dat/fourier_transform_ex_gauss_freq_1.0.txt' with lines,\
'dat/fourier_transform_ex_gauss_freq_2.0.txt' with lines
unset multiplot
set terminal pngcairo size 560,480 enhanced font 'Verdana,8'
set output 'img/fourier_series_pimp.png'
set multiplot layout 2,1 rowsfirst
#
set xlabel 't, c'
set ylabel 's(t)'
plot[-10:10] 'dat/fourier_transform_ex_gauss_time_0.5.txt' with lines,\
'dat/fourier_transform_ex_gauss_time_1.0.txt' with lines,\
'dat/fourier_transform_ex_gauss_time_2.0.txt' with lines
#
set xlabel 'w, рад/c'
set ylabel 'S(w)'
plot[-4*pi:4*pi] 'dat/fourier_transform_ex_gauss_freq_0.5.txt' with lines,\
'dat/fourier_transform_ex_gauss_freq_1.0.txt' with lines,\
'dat/fourier_transform_ex_gauss_freq_2.0.txt' with lines
unset multiplot

8
test/src/cheby2_ap_zp.c
Wyświetl plik

@ -23,11 +23,11 @@ int main()
printf("\nChebyshev type 2 zeros:\n");
for(k = 0; k < nz; k++)
printf("(%9.3f, %9.3f)\n", k, SCALE*RE(z[k]), SCALE*IM(z[k]));
printf("(%9.3f, %9.3f)\n", SCALE*RE(z[k]), SCALE*IM(z[k]));
printf("\nChebyshev type 2 poles:\n");
for(k = 0; k < np; k++)
printf("(%9.3f, %9.3f)\n", k, SCALE*RE(p[k]), SCALE*IM(p[k]));
printf("(%9.3f, %9.3f)\n", SCALE*RE(p[k]), SCALE*IM(p[k]));
double alpha[N], sigma[N], omega[N];
@ -67,11 +67,11 @@ int main()
printf("\nChebyshev type 2 zeros:\n");
for(k = 0; k < nz; k++)
printf("(%9.3f, %9.3f)\n", k, SCALE*RE(z[k]), SCALE*IM(z[k]));
printf("(%9.3f, %9.3f)\n", SCALE*RE(z[k]), SCALE*IM(z[k]));
printf("\nChebyshev type 2 poles:\n");
for(k = 0; k < np; k++)
printf("(%9.3f, %9.3f)\n", k, SCALE*RE(p[k]), SCALE*IM(p[k]));
printf("(%9.3f, %9.3f)\n", SCALE*RE(p[k]), SCALE*IM(p[k]));
for(m = 1; m < 12; m++)
{

42
test/src/fftn_create_test.c
Wyświetl plik

@ -1,42 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define N 7
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);
for(n = 0; n < N; n++)
{
printf("y[%3d] = %12.4f%12.4f\n", n, RE(y[n]), IM(y[n]));
}
fft_free(&pfft); // clear FFT object
dspl_free(handle); // free dspl handle
return 0;
}

70
test/src/fftn_performance.c
Wyświetl plik

@ -1,70 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include "dspl.h"
#define N 1048576
int main()
{
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
int n, k, m;
complex_t *x = NULL;
complex_t *y = NULL; // DFT
fft_t pfft; // FFT object
clock_t t0, t1;
double dt;
x = (complex_t*)malloc(N*sizeof(complex_t));
y = (complex_t*)malloc(N*sizeof(complex_t));
//clear fft object
memset(&pfft, 0, sizeof(fft_t));
for(n = 0; n < N; n++)
{
RE(x[n]) = (double)n;
IM(x[n]) = 0.0;
}
n = N;
k = 4;
while(n > 4)
{
printf("n= %16d ", n);
fftn_create(&pfft, n);
t0 = clock();
for(m = 0; m < k; m++)
fftn_krn(x, y, &pfft, n, 0);
t1 = clock();
dt = (1000.0 * (double) (t1 - t0)) / CLOCKS_PER_SEC / (double)m;
printf("%12.6f ms ", dt);
fft_free(&pfft); // clear FFT object
fft_create(&pfft, n);
t0 = clock();
for(m = 0; m < k; m++)
fft_cmplx(x, n, &pfft, y);
t1 = clock();
dt = (1000.0 * (double) (t1 - t0)) / CLOCKS_PER_SEC / (double)m;
printf("%12.6f ms \n", dt);
fft_free(&pfft); // clear FFT object
n/=2;
k = (int)(k*1.5);
}
dspl_free(handle); // free dspl handle
free(x);
free(y);
return 0;
}

149
test/src/filter_approx.c
Wyświetl plik

@ -3,126 +3,43 @@
#include <math.h>
#include "dspl.h"
#define N 8000
void r2hdb(double* r, double* h, double* hdb, double ep2, int n)
{
int k;
for(k = 0; k < n; k++)
{
r[k] *= r[k];
h[k] = 1.0/(1.0 + ep2 * r[k]);
hdb[k] = 10.0*log10(h[k]);
}
}
#define N 1000
int main(int argc, char* argv[])
{
double w[N]; // время (сек)
double r[N]; // входной сигнал
double h[N], hdb[N];
double ep2, Rp, Rs, es2, m;
int ord, k;
void* handle; // DSPL handle
handle = dspl_load(); // Load DSPL function
// заполняю массив частот в логарифмическом формате от 0.01 до 100
logspace(-2, 2, N, DSPL_PERIODIC, w);
ord = 4; // порядок фильтра
Rp = 0.3; // Неравномерность в полосе пропускания (дБ)
Rs = 15.0; // Уровень подавления в полосе заграждения (дБ)
// Параметры ep^2 и es^2
ep2 = pow(10.0, Rp*0.1) - 1.0;
es2 = pow(10.0, Rs*0.1) - 1.0;
// вывод на печать параметров ep^2 и es^2
printf("ep^2 = %.4f\n", ep2);
printf("es^2 = %.4f\n", es2);
//**********************************************************************
// Расчет F_N^2(w) и |H(jw)|^2 фильтра Баттерворта
//**********************************************************************
for(k = 0; k < N; k++)
{
r[k] = pow(w[k], (double)(ord));
}
r2hdb(r, h, hdb, ep2, N);
// сохранение в файлы результатов расчета фильтра Баттерворта
writetxt(w, r, N, "dat/butter_r.txt");
writetxt(w, h, N, "dat/butter_h.txt");
writetxt(w, hdb, N, "dat/butter_hdb.txt");
double w[N], r[N], h[N];
double ep2, Gp2;
int ord, k;
void* handle;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
//**********************************************************************
// Расчет F_N^2(w) и |H(jw)|^2 фильтра Чебышева 1-го рода
//**********************************************************************
cheby_poly1(w, N, ord, r);
r2hdb(r, h, hdb, ep2, N);
// сохранение в файлы результатов расчета фильтра Чебышева 1-го рода
writetxt(w, r, N, "dat/cheby1_r.txt");
writetxt(w, h, N, "dat/cheby1_h.txt");
writetxt(w, hdb, N, "dat/cheby1_hdb.txt");
//**********************************************************************
// Расчет F_N^2(w) и |H(jw)|^2 фильтра Чебышева 2-го рода
//**********************************************************************
for(k = 0; k < N; k++)
{
w[k] = 1.0 / w[k];
}
cheby_poly1(w, N, ord, r);
for(k = 0; k < N; k++)
{
r[k] =1.0 / (r[k] *r[k]);
h[k] = 1.0/(1.0 + es2*r[k]);
hdb[k] = 10.0*log10(h[k]);
}
logspace(-2, 2, N, DSPL_PERIODIC, w);
// сохранение в файлы результатов расчета фильтра Чебышева 1-го рода
writetxt(w, r, N, "dat/cheby2_r.txt");
writetxt(w, h, N, "dat/cheby2_h.txt");
writetxt(w, hdb, N, "dat/cheby2_hdb.txt");
//**********************************************************************
// Расчет F_N^2(w) и |H(jw)|^2 эллиптического фильтра
//**********************************************************************
ellip_modulareq(Rp, Rs, ord, &m); // пересчет эллиптического модуля
printf("modular m = %.3f\n", m); // вывод на печать
// расчет эллиптической рациональной функции
ellip_rat(w, N, ord, m, r);
r2hdb(r, h, hdb, ep2, N);
// сохранение в файлы результатов расчета эллиптического фильтра
writetxt(w, r, N, "dat/ellip_r.txt");
writetxt(w, h, N, "dat/ellip_h.txt");
writetxt(w, hdb, N, "dat/ellip_hdb.txt");
linspace(0, 2.5, N, DSPL_PERIODIC, w);
dspl_free(handle); // free dspl handle
return 0;
ord = 4;
Gp2 = 0.9;
ep2 = 1.0 / Gp2 -1;
for(k = 0; k < N; k++)
{
r[k] = pow(w[k], (double)(ord));
r[k] *= r[k];
h[k] = 6.0/(1.0 + ep2*r[k]);
w[k] *= 4;
}
writetxt(w, h, N, "dat/butter_approx.txt");
// remember to free the resource
dspl_free(handle);
return 0;
}

93
test/src/fourier_series_pimp_q.c 100755
Wyświetl plik

@ -0,0 +1,93 @@
#include <stdio.h>
#include <string.h>
//#include "common.h"
#include "dspl.h"
//#include "plot.h"
#define N 1000
#define T 4.0
#define A 2.0
#define M 41
int main(int argc, char* argv[])
{
double t[N]; // время (сек)
double s[N]; // входной сигнал
complex_t S[M]; // комплексный спектр периодического сигнала
double Smag[M]; // амплитудный спектр периодического сигнала
double w[M]; // частота (рад/c) дискретного спектра
double wc[N]; // частота (рад/с) огибающей спектра
double Sc[N]; // огибающая спектра
double tau; // длительность импульса
// скважность
double Q[3] = {5.0, 2.0, 1.25};
int q, m, n;
char fname[64];
void* handle;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
for(q = 0; q < 3; q++)
{
tau = T/Q[q];
// заполнение массива временных отсчетов
// на 4-x периодах повторения сигнала
// для отображения на осциллограмме
linspace(-T*2.0, T*2.0, N, DSPL_PERIODIC, t);
// 4 периода повторения п-импульса скважности Q[q]
signal_pimp(t, N, A, tau, 0.0, T, s);
// сохранение в текстовый файл
sprintf(fname, "dat/pimp_time_%.2lf.csv", Q[q]);
writetxt(t, s, N, fname);
// заполнение массива временных отсчетов
// на одном периоде повторения сигнала
linspace(-T/2.0, T/2.0, N, DSPL_PERIODIC, t);
// один период повторения п-импульса скважности Q[q]
signal_pimp(t, N, A, tau, 0.0, T, s);
// разложение в ряд Фурье
fourier_series_dec(t, s, N, T, M, w, S);
// Рассчет амплитудного спектра
for(m = 0; m < M; m++)
{
printf("S[%d] = %f %f\n", m, RE(S[m]), IM(S[m]));
Smag[m] = ABS(S[m]);
}
// Сохранение в файл амплитудного спетра для скважности Q[q]
sprintf(fname, "dat/pimp_freq_discrete_%.2lf.csv", Q[q]);
writetxt(w, Smag, M, fname);
// Вектор частот непрерывной огибаюхей вида sin(w/2*tau) / (w/2*T)
linspace(w[0], w[M-1], N, DSPL_SYMMETRIC, wc);
// Расчет огибающей
for(n = 0; n < N; n++)
Sc[n] = (wc[n] == 0.0) ? A/Q[q] : fabs( A * sin(0.5*wc[n]*tau) / (0.5*wc[n] * T));
// сохранение огибающей в файл для скважности Q[q]
sprintf(fname, "dat/pimp_freq_cont_%.2lf.csv", Q[q]);
writetxt(wc, Sc, N, fname);
}
// remember to free the resource
dspl_free(handle);
return 0;
}

4
test/src/fourier_series_pimp_spectrum.c
Wyświetl plik

@ -1,4 +1,5 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "dspl.h"
@ -17,7 +18,8 @@ int main(int argc, char* argv[])
double mag[K]; // амплитудный спектр
double phi[K]; // фазовый спектр
void* handle;
int k, n;
// test
int n;
handle = dspl_load();
if(!handle)

114
test/src/fourier_transform_ex.c 100755
Wyświetl plik

@ -0,0 +1,114 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "dspl.h"
#define N 2000
#define A 1
int main(int argc, char* argv[])
{
double t[N], w[N], s[N], S[N], PHI[N], sigma;
void* handle;
int k;
handle = dspl_load();
if(!handle)
{
printf("cannot to load libdspl!\n");
return 0;
}
linspace(-10, 10, N, DSPL_PERIODIC, t);
linspace(-4*M_PI, 4*M_PI, N, DSPL_PERIODIC, w);
sigma = 0.5;
for(k = 0; k < N; k++)
{
s[k] = A * exp(-t[k]*t[k]/(sigma*sigma));
S[k] = A * sigma * sqrt(M_PI) * exp(-w[k]*w[k]*sigma*sigma*0.25);
}
writetxt(t, s, N, "dat/fourier_transform_ex_gauss_time_0.5.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_gauss_freq_0.5.txt");
for(k = 0; k < N; k++)
{
s[k] = A * exp(-fabs(t[k]) * sigma);
S[k] = 2 * A * sigma / (sigma*sigma + w[k]*w[k]);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp2_time_0.5.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp2_freq_0.5.txt");
for(k = 0; k < N; k++)
{
s[k] = t[k] >=0 ? A * exp(-t[k] * sigma) : 0.0;
S[k] = A / sqrt(sigma*sigma + w[k]*w[k]);
PHI[k] = -atan(w[k] / sigma);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp1_time_0.5.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp_freq_mag_0.5.txt");
writetxt(w, PHI, N, "dat/fourier_transform_ex_exp_freq_phi_0.5.txt");
sigma = 1;
for(k = 0; k < N; k++)
{
s[k] = A * exp(-t[k]*t[k]/(sigma*sigma));
S[k] = A * sigma * sqrt(M_PI) * exp(-w[k]*w[k]*sigma*sigma*0.25);
}
writetxt(t, s, N, "dat/fourier_transform_ex_gauss_time_1.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_gauss_freq_1.0.txt");
for(k = 0; k < N; k++)
{
s[k] = A * exp(-fabs(t[k]) * sigma);
S[k] = 2 * A * sigma / (sigma*sigma + w[k]*w[k]);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp2_time_1.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp2_freq_1.0.txt");
for(k = 0; k < N; k++)
{
s[k] = t[k] >=0 ? A * exp(-t[k] * sigma) : 0.0;
S[k] = A / sqrt(sigma*sigma + w[k]*w[k]);
PHI[k] = -atan(w[k] / sigma);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp1_time_1.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp_freq_mag_1.0.txt");
writetxt(w, PHI, N, "dat/fourier_transform_ex_exp_freq_phi_1.0.txt");
sigma = 2;
for(k = 0; k < N; k++)
{
s[k] = A * exp(-t[k]*t[k]/(sigma*sigma));
S[k] = A * sigma * sqrt(M_PI) * exp(-w[k]*w[k]*sigma*sigma*0.25);
}
writetxt(t, s, N, "dat/fourier_transform_ex_gauss_time_2.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_gauss_freq_2.0.txt");
for(k = 0; k < N; k++)
{
s[k] = A * exp(-fabs(t[k]) * sigma);
S[k] = 2 * A * sigma / (sigma*sigma + w[k]*w[k]);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp2_time_2.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp2_freq_2.0.txt");
for(k = 0; k < N; k++)
{
s[k] = t[k] >=0 ? A * exp(-t[k] * sigma) : 0.0;
S[k] = A / sqrt(sigma*sigma + w[k]*w[k]);
PHI[k] = -atan(w[k] / sigma);
}
writetxt(t, s, N, "dat/fourier_transform_ex_exp1_time_2.0.txt");
writetxt(w, S, N, "dat/fourier_transform_ex_exp_freq_mag_2.0.txt");
writetxt(w, PHI, N, "dat/fourier_transform_ex_exp_freq_phi_2.0.txt");
// remember to free the resource
dspl_free(handle);
return system("gnuplot -p gnuplot/fourier_transform_ex_gauss.plt");;
}

3
test/src/fourier_transform_period_incr.c
Wyświetl plik

@ -1,5 +1,6 @@
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include "dspl.h"
@ -16,7 +17,7 @@ int main(int argc, char* argv[])
double w[K], mag[K];
double T;
void* handle;
int k, n;
int n;

2
test/src/laplace_tf_planes.c
Wyświetl plik

@ -15,7 +15,7 @@ int main(int argc, char* argv[])
double w[N], sigma[N], hr[N], hi[N];
void* handle;
double b[ORD+1], a[ORD+1];
complex_t p[ORD], z[ORD];
//complex_t p[ORD], z[ORD];
complex_t hs[N];
int k;

23
verif/bin/octave/functions/dspl_readbin.m 100755
Wyświetl plik

@ -0,0 +1,23 @@
function [dat, n, m] = dspl_readbin(fn)
fid = fopen(fn);
if(~fid)
error('cannot to open file');
end
type = fread(fid, 1, 'int32');
n = fread(fid, 1, 'int32');
m = fread(fid, 1, 'int32');
if(type==0)
dat = fread(fid, [n*m, 1], 'double');
dat = reshape(dat, n, m);
end
if(type==1)
dat = fread(fid, [n*m*2, 1], 'double');
dat = dat(1:2:end) + 1i * dat(2:2:end);
dat = reshape(dat, n, m);
end
fclose(fid);
end

35
verif/bin/octave/functions/dspl_writebin.m 100755
Wyświetl plik

@ -0,0 +1,35 @@
function dspl_writebin(dat, fn)
fid = fopen(fn, 'w');
if(~fid)
error('cannot to open file');
end
if(isreal (dat))
type = 0;
else
type = 1;
end
n = size(dat,1);
m = size(dat,2);
fwrite (fid, type, 'int32');
fwrite (fid, n, 'int32');
fwrite (fid, m, 'int32');
if(type==0)
dat = reshape(dat, n*m, 1);
fwrite (fid, dat, 'double');
end
if(type==1)
y = zeros(2*n*m, 1);
dat = reshape(dat, n*m, 1);
y(1:2:end) = real(dat);
y(2:2:end) = imag(dat);
fwrite (fid, y, 'double');
end
fclose(fid);
end

6
verif/bin/octave/readbin_verif.m 100755
Wyświetl plik

@ -0,0 +1,6 @@
clear all; close all; clc;
addpath('octave/functions');
x = dspl_readbin('dat/in.dat');
dspl_writebin(x, 'dat/out.dat');

23
verif/src/readbin_verif.c
Wyświetl plik

@ -20,17 +20,28 @@ int main()
randn(xr, N, 0, 1.0);
randn((double*)xc, 2*N, 0, 1.0);
writebin(xr, N, DAT_DOUBLE, "dat/in_real.dat");
writebin(xc, N, DAT_COMPLEX, "dat/in_cmplx.dat");
readbin("dat/in_real.dat", (void**)&yr, NULL, NULL);
readbin("dat/in_cmplx.dat", (void**)&yc, NULL, NULL);
writebin(xr, N, DAT_DOUBLE, "dat/in.dat");
system("octave octave/readbin_verif.m");
readbin("dat/out.dat", (void**)&yr, NULL, NULL);
vr = verif(xr, yr, N, 1E-12, &err);
printf("readbin real verification error: %12.4e\n", err);
if(vr == DSPL_VERIF_SUCCESS)
printf("readbin real verification OK: %12.4e\n", err);
else
printf("readbin real verification ERROR: %12.4e\n", err);
writebin(xc, N, DAT_COMPLEX, "dat/in.dat");
system("octave octave/readbin_verif.m");
readbin("dat/out.dat", (void**)&yc, NULL, NULL);
vr = verif_cmplx(xc, yc, N, 1E-12, &err);
printf("readbin cmplx verification error: %12.4e\n", err);
if(vr == DSPL_VERIF_SUCCESS)
printf("readbin cmplx verification OK: %12.4e\n", err);
else
printf("readbin cmplx verification ERROR: %12.4e\n", err);
dspl_free(handle); // free dspl handle