/*******************************************************************************
* @brief One step 1D 5/3 discrete wavelet transform.
* @param in The input data array.
* @param out The output data array.
* @param w The length
*******************************************************************************/
void dwt_53_1d( int16_t *in, int16_t *out, const uint32_t w )
{
int wt = w - 2, i, j;
int16_t *l, *h;
l = out; h = &out[(w>>1) + (w&1)];
h[0] = in[1] - ((in[0] + in[2])>>1);
l[0] = in[0] + (h[0]>>1);
for(i=2,j=1; i < wt; i+=2,j++)
{
h[j] = in[i+1] - ((in[i] + in[i+2])>>1);
l[j] = in[i] + ((h[j-1] + h[j])>>2);
}
if(w&1)
{
l[j] = in[i] + (h[j-1]>>1);
}
else
{
h[j] = in[i+1] - in[i];
l[j] = in[i] + ((h[j-1] + h[j])>>2);
}
}
/*******************************************************************************
* @brief >Compute fast Fourier transform
* @param >signal buffer for transform
* @param >output data
* @param >buffer length
* @param >1 - forward transform, -1 backward transform
* @note >buffer length must be equal 64,128,256 etc.
* @verbm >re im и к-во точек исходного массива. потом 1 - прямое, -1 обратное преобразование
* есть ограничение - массив должен быть кратен 2 - лучше 64,128,256 и т.д
* re это массив который вводится (реальная часть), im можно обнулить просто (вещественная часть)
*******************************************************************************/
void fft( float *re, float *im, unsigned long n, int isign )
{
unsigned long i, j, k, l, le, le1, ip, n2;
float wpr, wpi, wr, wi, wtr, wti;
n2 = n>>1;
j = 1;
for(i=0; i<n-1; i++)
{
if(i<j)
{
wtr = re[j-1];
wti = im[j-1];
re[j-1] = re[i];
im[j-1] = im[i];
re[i] = wtr;
im[i] = wti;
}
k = n2;
while(k<j){
j -= k;
k >>= 1;
}
j += k;
}
l=1;
k=n;
while(k>>=1)
{
le1 = (le=1<<l++) >> 1;
wtr = PI / (float)le1;
wpr = cos(wtr); wpi = -isign*sin(wtr);
wr = 1.0; wi = 0.0;
for(j=0; j<le1; j++)
{
for(i=j; i<n; i+=le)
{
ip = i + le1;
wtr = wr*re[ip] - wi*im[ip];
wti = wi*re[ip] + wr*im[ip];
re[ip] = re[i] - wtr;
im[ip] = im[i] - wti;
re[i] = re[i] + wtr;
im[i] = im[i] + wti;
}
wr = (wtr=wr)*wpr - wi*wpi;
wi = wi*wpr + wtr*wpi;
}
}
}
Средний
