35 changed files with 5121 additions and 29 deletions
@ -0,0 +1,169 @@ |
|||||||
|
// speex_aec.c — C-обёртка над SpeexDSP mdf.c (акустическое эхоподавление).
|
||||||
|
//
|
||||||
|
// API см. в speex_aec.h. Ключевые детали реализации:
|
||||||
|
// - собственная линия задержки рендера на depth кадров (SpeexDSP даёт только
|
||||||
|
// фиксированные 2 кадра — PLAYBACK_DELAY в mdf.c);
|
||||||
|
// - синхронный speex_echo_cancellation(rec, play_delayed, out) с выравниванием
|
||||||
|
// рендер↔захват по линии задержки;
|
||||||
|
// - дрейф независимых потоков: переполнение линии → отбрасываем старый кадр,
|
||||||
|
// недозаполнение (старт/underflow) → passthrough без канселлера.
|
||||||
|
//
|
||||||
|
// Захват и рендер приходят чанками разного размера (480/960), поэтому обе
|
||||||
|
// стороны накапливаются до полного кадра (frame_samples) внутри обёртки.
|
||||||
|
|
||||||
|
#include "speex_aec.h" |
||||||
|
#include "debug_config.h" |
||||||
|
#include "mem.h" |
||||||
|
|
||||||
|
#include <string.h> |
||||||
|
|
||||||
|
#include <speex/speex_echo.h> |
||||||
|
|
||||||
|
#define AEC_ID "speex_aec" |
||||||
|
|
||||||
|
struct speex_aec { |
||||||
|
SpeexEchoState* st; |
||||||
|
int frame_samples; /* сэмплов в кадре (960 @48кГц) */ |
||||||
|
int depth; /* глубина линии задержки в кадрах (>=1) */ |
||||||
|
int16_t* line; /* кольцо depth*frame_samples */ |
||||||
|
int head; /* индекс самого старого кадра в line */ |
||||||
|
int count; /* кадров в линии (0..depth) */ |
||||||
|
int16_t* play_acc; /* аккумулятор рендера до кадра */ |
||||||
|
int play_len; |
||||||
|
int16_t* cap_acc; /* аккумулятор захвата до кадра */ |
||||||
|
int cap_len; |
||||||
|
uint32_t overruns; /* отброшено рендер-кадров (дрейф: рендер быстрее) */ |
||||||
|
uint32_t underruns; /* passthrough-кадров (дрейф: рендер медленнее/старт) */ |
||||||
|
}; |
||||||
|
|
||||||
|
/* Положить полный рендер-кадр в линию задержки (вытесняя старейший при переполнении). */ |
||||||
|
static void aec_enqueue(speex_aec_t* a, const int16_t* frame) { |
||||||
|
if (a->count == a->depth) { |
||||||
|
a->head = (a->head + 1) % a->depth; |
||||||
|
a->count--; |
||||||
|
a->overruns++; |
||||||
|
if (a->overruns == 1 || (a->overruns & 0xff) == 0) { |
||||||
|
DEBUG_WARN(DEBUG_CATEGORY_AEC, |
||||||
|
"%s: delay line overflow (render faster than capture) overruns=%u", AEC_ID, a->overruns); |
||||||
|
} |
||||||
|
} |
||||||
|
memcpy(a->line + (size_t)((a->head + a->count) % a->depth) * a->frame_samples, |
||||||
|
frame, (size_t)a->frame_samples * sizeof(int16_t)); |
||||||
|
a->count++; |
||||||
|
} |
||||||
|
|
||||||
|
speex_aec_t* speex_aec_create(int sample_rate, int frame_samples, int filter_samples, int delay_frames) { |
||||||
|
speex_aec_t* a; |
||||||
|
int rate; |
||||||
|
|
||||||
|
if (sample_rate <= 0 || frame_samples <= 0 || filter_samples < frame_samples) { |
||||||
|
DEBUG_ERROR(DEBUG_CATEGORY_AEC, "%s: bad args rate=%d frame=%d filter=%d delay=%d", |
||||||
|
AEC_ID, sample_rate, frame_samples, filter_samples, delay_frames); |
||||||
|
return NULL; |
||||||
|
} |
||||||
|
|
||||||
|
a = (speex_aec_t*)u_calloc(1, sizeof(*a)); |
||||||
|
if (!a) { |
||||||
|
DEBUG_ERROR(DEBUG_CATEGORY_AEC, "%s: OOM for state", AEC_ID); |
||||||
|
return NULL; |
||||||
|
} |
||||||
|
|
||||||
|
a->st = speex_echo_state_init(frame_samples, filter_samples); |
||||||
|
rate = sample_rate; |
||||||
|
if (speex_echo_ctl(a->st, SPEEX_ECHO_SET_SAMPLING_RATE, &rate) != 0) { |
||||||
|
DEBUG_ERROR(DEBUG_CATEGORY_AEC, "%s: SET_SAMPLING_RATE(%d) failed", AEC_ID, sample_rate); |
||||||
|
speex_echo_state_destroy(a->st); |
||||||
|
u_free(a); |
||||||
|
return NULL; |
||||||
|
} |
||||||
|
|
||||||
|
a->frame_samples = frame_samples; |
||||||
|
a->depth = delay_frames > 0 ? delay_frames : 1; |
||||||
|
a->line = (int16_t*)u_calloc((uint32_t)(a->depth * frame_samples), sizeof(int16_t)); |
||||||
|
a->play_acc = (int16_t*)u_calloc((uint32_t)frame_samples, sizeof(int16_t)); |
||||||
|
a->cap_acc = (int16_t*)u_calloc((uint32_t)frame_samples, sizeof(int16_t)); |
||||||
|
if (!a->line || !a->play_acc || !a->cap_acc) { |
||||||
|
DEBUG_ERROR(DEBUG_CATEGORY_AEC, "%s: OOM for buffers", AEC_ID); |
||||||
|
if (a->line) u_free(a->line); |
||||||
|
if (a->play_acc) u_free(a->play_acc); |
||||||
|
if (a->cap_acc) u_free(a->cap_acc); |
||||||
|
speex_echo_state_destroy(a->st); |
||||||
|
u_free(a); |
||||||
|
return NULL; |
||||||
|
} |
||||||
|
|
||||||
|
DEBUG_INFO(DEBUG_CATEGORY_AEC, "%s: created rate=%d frame=%d filter=%d delay=%d frames", |
||||||
|
AEC_ID, sample_rate, frame_samples, filter_samples, a->depth); |
||||||
|
return a; |
||||||
|
} |
||||||
|
|
||||||
|
void speex_aec_destroy(speex_aec_t* a) { |
||||||
|
if (!a) return; |
||||||
|
speex_echo_state_destroy(a->st); |
||||||
|
u_free(a->line); |
||||||
|
u_free(a->play_acc); |
||||||
|
u_free(a->cap_acc); |
||||||
|
u_free(a); |
||||||
|
} |
||||||
|
|
||||||
|
void speex_aec_reset(speex_aec_t* a) { |
||||||
|
if (!a) return; |
||||||
|
speex_echo_state_reset(a->st); |
||||||
|
a->head = 0; |
||||||
|
a->count = 0; |
||||||
|
a->play_len = 0; |
||||||
|
a->cap_len = 0; |
||||||
|
DEBUG_INFO(DEBUG_CATEGORY_AEC, "%s: reset (filter + delay line)", AEC_ID); |
||||||
|
} |
||||||
|
|
||||||
|
void speex_aec_feed_playback(speex_aec_t* a, const int16_t* pcm, int count) { |
||||||
|
int fs, take; |
||||||
|
|
||||||
|
if (!a || !pcm || count <= 0) return; |
||||||
|
fs = a->frame_samples; |
||||||
|
|
||||||
|
take = count; |
||||||
|
if (a->play_len + count > fs) take = fs - a->play_len; |
||||||
|
memcpy(a->play_acc + a->play_len, pcm, (size_t)take * sizeof(int16_t)); |
||||||
|
a->play_len += take; |
||||||
|
if (a->play_len < fs) return; |
||||||
|
|
||||||
|
aec_enqueue(a, a->play_acc); |
||||||
|
a->play_len = 0; |
||||||
|
} |
||||||
|
|
||||||
|
int speex_aec_process_capture(speex_aec_t* a, const int16_t* pcm, int count, int16_t* out) { |
||||||
|
int fs, take; |
||||||
|
|
||||||
|
if (!a || !pcm || !out || count <= 0) return 0; |
||||||
|
fs = a->frame_samples; |
||||||
|
|
||||||
|
take = count; |
||||||
|
if (a->cap_len + count > fs) take = fs - a->cap_len; |
||||||
|
memcpy(a->cap_acc + a->cap_len, pcm, (size_t)take * sizeof(int16_t)); |
||||||
|
a->cap_len += take; |
||||||
|
if (a->cap_len < fs) return 0; |
||||||
|
|
||||||
|
if (a->count == a->depth) { |
||||||
|
const int16_t* ref = a->line + (size_t)a->head * fs; |
||||||
|
speex_echo_cancellation(a->st, a->cap_acc, ref, out); |
||||||
|
a->head = (a->head + 1) % a->depth; |
||||||
|
a->count--; |
||||||
|
} else { |
||||||
|
/* линия не наполнена (старт или рендер отстаёт): passthrough без канселлера */ |
||||||
|
if (a->underruns == 0 || (a->underruns & 0xff) == 0) { |
||||||
|
DEBUG_WARN(DEBUG_CATEGORY_AEC, |
||||||
|
"%s: delay line underflow (render slower than capture) underruns=%u fill=%d/%d", |
||||||
|
AEC_ID, a->underruns, a->count, a->depth); |
||||||
|
} |
||||||
|
memcpy(out, a->cap_acc, (size_t)fs * sizeof(int16_t)); |
||||||
|
a->underruns++; |
||||||
|
} |
||||||
|
|
||||||
|
a->cap_len = 0; |
||||||
|
return fs; |
||||||
|
} |
||||||
|
|
||||||
|
int speex_aec_delay_fill(const speex_aec_t* a) { |
||||||
|
return a ? a->count : 0; |
||||||
|
} |
||||||
@ -0,0 +1,74 @@ |
|||||||
|
/*
|
||||||
|
* speex_aec.h — акустическое эхоподавление (AEC), C-обёртка над SpeexDSP mdf.c. |
||||||
|
* |
||||||
|
* Подавляет эхо дальнего конца (то, что играем в динамик) в сигнале микрофона |
||||||
|
* перед кодированием. Работает на int16 PCM, один канал, частота 48000 (можно |
||||||
|
* 8000/16000/32000/48000), кадр 20 мс (960 сэмплов @48 кГц). |
||||||
|
* |
||||||
|
* Модель использования (duplex-контур звонка): |
||||||
|
* - рендер (far-end, то что пошло в динамик) → speex_aec_feed_playback(); |
||||||
|
* - захват (near-end, микрофон) → speex_aec_process_capture(). |
||||||
|
* |
||||||
|
* Обёртка держит собственную линию задержки рендера на `delay_frames` кадров и |
||||||
|
* зовёт синхронный speex_echo_cancellation(rec, play_delayed, out). Это делает |
||||||
|
* выравнивание рендер↔захват предсказуемым и настраиваемым (встроенный буфер |
||||||
|
* SpeexDSP фиксирован в 2 кадра и для Android-задержки не годится). |
||||||
|
* |
||||||
|
* Дрейф двух независимых потоков (capture/play на Android) компенсируется |
||||||
|
* ограниченной глубиной линии: переполнение → отбрасываем старый рендер-кадр, |
||||||
|
* недозаполнение → passthrough без канселлера (с подробным логом категории "aec"). |
||||||
|
* |
||||||
|
* Зависимости: SpeexDSP (mdf.c/fftwrap.c/kiss_fft*) из lib/speexdsp, флаги |
||||||
|
* FLOATING_POINT + USE_KISS_FFT. Собирается всегда (без внешних зависимостей). |
||||||
|
* |
||||||
|
* Лицензия SpeexDSP: 3-clause BSD (Xiph) — см. lib/speexdsp/COPYING. |
||||||
|
*/ |
||||||
|
#ifndef SPEEX_AEC_H |
||||||
|
#define SPEEX_AEC_H |
||||||
|
|
||||||
|
#include <stdint.h> |
||||||
|
|
||||||
|
#ifdef __cplusplus |
||||||
|
extern "C" { |
||||||
|
#endif |
||||||
|
|
||||||
|
typedef struct speex_aec speex_aec_t; |
||||||
|
|
||||||
|
/**
|
||||||
|
* Создать эхоканселлер. |
||||||
|
* sample_rate — 48000 (также 8000/16000/32000); |
||||||
|
* frame_samples — сэмплов в кадре (960 = 20 мс @48 кГц); |
||||||
|
* filter_samples— длина эхо-хвоста в сэмплах (14400 = 300 мс, кратно кадру); |
||||||
|
* delay_frames — задержка рендер→захват в кадрах (>=0; 0 = без линии задержки). |
||||||
|
* Возвращает NULL при ошибке (лог категории "aec"). |
||||||
|
*/ |
||||||
|
speex_aec_t* speex_aec_create(int sample_rate, int frame_samples, int filter_samples, int delay_frames); |
||||||
|
|
||||||
|
/* Освободить канселлер. NULL безопасен. */ |
||||||
|
void speex_aec_destroy(speex_aec_t* aec); |
||||||
|
|
||||||
|
/* Сбросить адаптивный фильтр и линию задержки (смена устройства/роута). */ |
||||||
|
void speex_aec_reset(speex_aec_t* aec); |
||||||
|
|
||||||
|
/**
|
||||||
|
* Подать рендер (дальний конец). Накопление до кадра внутри; полный кадр |
||||||
|
* кладётся в линию задержки. count может быть любым (обычно 480 или 960). |
||||||
|
*/ |
||||||
|
void speex_aec_feed_playback(speex_aec_t* aec, const int16_t* pcm, int count); |
||||||
|
|
||||||
|
/**
|
||||||
|
* Обработать захват (микрофон). Накопление до кадра внутри. Когда кадр собран — |
||||||
|
* подавляет эхо (с учётом линии задержки) и пишет результат в out. |
||||||
|
* Возвращает число сэмплов, записанных в out (frame_samples, если кадр готов; |
||||||
|
* 0 — кадр ещё накапливается). out не должен алиаситься с pcm. |
||||||
|
*/ |
||||||
|
int speex_aec_process_capture(speex_aec_t* aec, const int16_t* pcm, int count, int16_t* out); |
||||||
|
|
||||||
|
/* Актуальная глубина линии задержки (кадров) — для диагностики дрейфа. */ |
||||||
|
int speex_aec_delay_fill(const speex_aec_t* aec); |
||||||
|
|
||||||
|
#ifdef __cplusplus |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
#endif /* SPEEX_AEC_H */ |
||||||
@ -0,0 +1,35 @@ |
|||||||
|
Copyright 2002-2008 Xiph.org Foundation |
||||||
|
Copyright 2002-2008 Jean-Marc Valin |
||||||
|
Copyright 2005-2007 Analog Devices Inc. |
||||||
|
Copyright 2005-2008 Commonwealth Scientific and Industrial Research |
||||||
|
Organisation (CSIRO) |
||||||
|
Copyright 1993, 2002, 2006 David Rowe |
||||||
|
Copyright 2003 EpicGames |
||||||
|
Copyright 1992-1994 Jutta Degener, Carsten Bormann |
||||||
|
|
||||||
|
Redistribution and use in source and binary forms, with or without |
||||||
|
modification, are permitted provided that the following conditions |
||||||
|
are met: |
||||||
|
|
||||||
|
- Redistributions of source code must retain the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer. |
||||||
|
|
||||||
|
- Redistributions in binary form must reproduce the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer in the |
||||||
|
documentation and/or other materials provided with the distribution. |
||||||
|
|
||||||
|
- Neither the name of the Xiph.org Foundation nor the names of its |
||||||
|
contributors may be used to endorse or promote products derived from |
||||||
|
this software without specific prior written permission. |
||||||
|
|
||||||
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
||||||
|
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
||||||
|
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
||||||
|
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR |
||||||
|
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
||||||
|
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
||||||
|
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
||||||
|
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
||||||
|
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
||||||
|
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
||||||
|
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
||||||
@ -0,0 +1,160 @@ |
|||||||
|
/*
|
||||||
|
Copyright (c) 2003-2004, Mark Borgerding |
||||||
|
|
||||||
|
All rights reserved. |
||||||
|
|
||||||
|
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: |
||||||
|
|
||||||
|
* Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. |
||||||
|
* Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. |
||||||
|
* Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission. |
||||||
|
|
||||||
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
||||||
|
*/ |
||||||
|
|
||||||
|
#define MIN(a,b) ((a)<(b) ? (a):(b)) |
||||||
|
#define MAX(a,b) ((a)>(b) ? (a):(b)) |
||||||
|
|
||||||
|
/* kiss_fft.h
|
||||||
|
defines kiss_fft_scalar as either short or a float type |
||||||
|
and defines |
||||||
|
typedef struct { kiss_fft_scalar r; kiss_fft_scalar i; }kiss_fft_cpx; */ |
||||||
|
#include "kiss_fft.h" |
||||||
|
#include "math_approx.h" |
||||||
|
|
||||||
|
#define MAXFACTORS 32 |
||||||
|
/* e.g. an fft of length 128 has 4 factors
|
||||||
|
as far as kissfft is concerned |
||||||
|
4*4*4*2 |
||||||
|
*/ |
||||||
|
|
||||||
|
struct kiss_fft_state{ |
||||||
|
int nfft; |
||||||
|
int inverse; |
||||||
|
int factors[2*MAXFACTORS]; |
||||||
|
kiss_fft_cpx twiddles[1]; |
||||||
|
}; |
||||||
|
|
||||||
|
/*
|
||||||
|
Explanation of macros dealing with complex math: |
||||||
|
|
||||||
|
C_MUL(m,a,b) : m = a*b |
||||||
|
C_FIXDIV( c , div ) : if a fixed point impl., c /= div. noop otherwise |
||||||
|
C_SUB( res, a,b) : res = a - b |
||||||
|
C_SUBFROM( res , a) : res -= a |
||||||
|
C_ADDTO( res , a) : res += a |
||||||
|
* */ |
||||||
|
#ifdef FIXED_POINT |
||||||
|
#include "arch.h" |
||||||
|
# define FRACBITS 15 |
||||||
|
# define SAMPPROD spx_int32_t |
||||||
|
#define SAMP_MAX 32767 |
||||||
|
|
||||||
|
#define SAMP_MIN -SAMP_MAX |
||||||
|
|
||||||
|
#if defined(CHECK_OVERFLOW) |
||||||
|
# define CHECK_OVERFLOW_OP(a,op,b) \ |
||||||
|
if ( (SAMPPROD)(a) op (SAMPPROD)(b) > SAMP_MAX || (SAMPPROD)(a) op (SAMPPROD)(b) < SAMP_MIN ) { \
|
||||||
|
fprintf(stderr,"WARNING:overflow @ " __FILE__ "(%d): (%d " #op" %d) = %ld\n",__LINE__,(a),(b),(SAMPPROD)(a) op (SAMPPROD)(b) ); } |
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
# define smul(a,b) ( (SAMPPROD)(a)*(b) ) |
||||||
|
# define sround( x ) (kiss_fft_scalar)( ( (x) + (1<<(FRACBITS-1)) ) >> FRACBITS ) |
||||||
|
|
||||||
|
# define S_MUL(a,b) sround( smul(a,b) ) |
||||||
|
|
||||||
|
# define C_MUL(m,a,b) \ |
||||||
|
do{ (m).r = sround( smul((a).r,(b).r) - smul((a).i,(b).i) ); \
|
||||||
|
(m).i = sround( smul((a).r,(b).i) + smul((a).i,(b).r) ); }while(0) |
||||||
|
|
||||||
|
# define C_MUL4(m,a,b) \ |
||||||
|
do{ (m).r = PSHR32( smul((a).r,(b).r) - smul((a).i,(b).i),17 ); \
|
||||||
|
(m).i = PSHR32( smul((a).r,(b).i) + smul((a).i,(b).r),17 ); }while(0) |
||||||
|
|
||||||
|
# define DIVSCALAR(x,k) \ |
||||||
|
(x) = sround( smul( x, SAMP_MAX/k ) ) |
||||||
|
|
||||||
|
# define C_FIXDIV(c,div) \ |
||||||
|
do { DIVSCALAR( (c).r , div); \
|
||||||
|
DIVSCALAR( (c).i , div); }while (0) |
||||||
|
|
||||||
|
# define C_MULBYSCALAR( c, s ) \ |
||||||
|
do{ (c).r = sround( smul( (c).r , s ) ) ;\
|
||||||
|
(c).i = sround( smul( (c).i , s ) ) ; }while(0) |
||||||
|
|
||||||
|
#else /* not FIXED_POINT*/ |
||||||
|
|
||||||
|
# define S_MUL(a,b) ( (a)*(b) ) |
||||||
|
#define C_MUL(m,a,b) \ |
||||||
|
do{ (m).r = (a).r*(b).r - (a).i*(b).i;\
|
||||||
|
(m).i = (a).r*(b).i + (a).i*(b).r; }while(0) |
||||||
|
|
||||||
|
#define C_MUL4(m,a,b) C_MUL(m,a,b) |
||||||
|
|
||||||
|
# define C_FIXDIV(c,div) /* NOOP */ |
||||||
|
# define C_MULBYSCALAR( c, s ) \ |
||||||
|
do{ (c).r *= (s);\
|
||||||
|
(c).i *= (s); }while(0) |
||||||
|
#endif |
||||||
|
|
||||||
|
#ifndef CHECK_OVERFLOW_OP |
||||||
|
# define CHECK_OVERFLOW_OP(a,op,b) /* noop */ |
||||||
|
#endif |
||||||
|
|
||||||
|
#define C_ADD( res, a,b)\ |
||||||
|
do { \
|
||||||
|
CHECK_OVERFLOW_OP((a).r,+,(b).r)\
|
||||||
|
CHECK_OVERFLOW_OP((a).i,+,(b).i)\
|
||||||
|
(res).r=(a).r+(b).r; (res).i=(a).i+(b).i; \
|
||||||
|
}while(0) |
||||||
|
#define C_SUB( res, a,b)\ |
||||||
|
do { \
|
||||||
|
CHECK_OVERFLOW_OP((a).r,-,(b).r)\
|
||||||
|
CHECK_OVERFLOW_OP((a).i,-,(b).i)\
|
||||||
|
(res).r=(a).r-(b).r; (res).i=(a).i-(b).i; \
|
||||||
|
}while(0) |
||||||
|
#define C_ADDTO( res , a)\ |
||||||
|
do { \
|
||||||
|
CHECK_OVERFLOW_OP((res).r,+,(a).r)\
|
||||||
|
CHECK_OVERFLOW_OP((res).i,+,(a).i)\
|
||||||
|
(res).r += (a).r; (res).i += (a).i;\
|
||||||
|
}while(0) |
||||||
|
|
||||||
|
#define C_SUBFROM( res , a)\ |
||||||
|
do {\
|
||||||
|
CHECK_OVERFLOW_OP((res).r,-,(a).r)\
|
||||||
|
CHECK_OVERFLOW_OP((res).i,-,(a).i)\
|
||||||
|
(res).r -= (a).r; (res).i -= (a).i; \
|
||||||
|
}while(0) |
||||||
|
|
||||||
|
|
||||||
|
#ifdef FIXED_POINT |
||||||
|
# define KISS_FFT_COS(phase) floor(MIN(32767,MAX(-32767,.5+32768 * cos (phase)))) |
||||||
|
# define KISS_FFT_SIN(phase) floor(MIN(32767,MAX(-32767,.5+32768 * sin (phase)))) |
||||||
|
# define HALF_OF(x) ((x)>>1) |
||||||
|
#elif defined(USE_SIMD) |
||||||
|
# define KISS_FFT_COS(phase) _mm_set1_ps( cos(phase) ) |
||||||
|
# define KISS_FFT_SIN(phase) _mm_set1_ps( sin(phase) ) |
||||||
|
# define HALF_OF(x) ((x)*_mm_set1_ps(.5)) |
||||||
|
#else |
||||||
|
# define KISS_FFT_COS(phase) (kiss_fft_scalar) cos(phase) |
||||||
|
# define KISS_FFT_SIN(phase) (kiss_fft_scalar) sin(phase) |
||||||
|
# define HALF_OF(x) ((x)*.5) |
||||||
|
#endif |
||||||
|
|
||||||
|
#define kf_cexp(x,phase) \ |
||||||
|
do{ \
|
||||||
|
(x)->r = KISS_FFT_COS(phase);\
|
||||||
|
(x)->i = KISS_FFT_SIN(phase);\
|
||||||
|
}while(0) |
||||||
|
#define kf_cexp2(x,phase) \ |
||||||
|
do{ \
|
||||||
|
(x)->r = spx_cos_norm((phase));\
|
||||||
|
(x)->i = spx_cos_norm((phase)-32768);\
|
||||||
|
}while(0) |
||||||
|
|
||||||
|
|
||||||
|
/* a debugging function */ |
||||||
|
#define pcpx(c)\ |
||||||
|
fprintf(stderr,"%g + %gi\n",(double)((c)->r),(double)((c)->i) ) |
||||||
@ -0,0 +1,232 @@ |
|||||||
|
/* Copyright (C) 2003 Jean-Marc Valin */ |
||||||
|
/**
|
||||||
|
@file arch.h |
||||||
|
@brief Various architecture definitions Speex |
||||||
|
*/ |
||||||
|
/*
|
||||||
|
Redistribution and use in source and binary forms, with or without |
||||||
|
modification, are permitted provided that the following conditions |
||||||
|
are met: |
||||||
|
|
||||||
|
- Redistributions of source code must retain the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer. |
||||||
|
|
||||||
|
- Redistributions in binary form must reproduce the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer in the |
||||||
|
documentation and/or other materials provided with the distribution. |
||||||
|
|
||||||
|
- Neither the name of the Xiph.org Foundation nor the names of its |
||||||
|
contributors may be used to endorse or promote products derived from |
||||||
|
this software without specific prior written permission. |
||||||
|
|
||||||
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
||||||
|
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
||||||
|
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
||||||
|
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR |
||||||
|
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
||||||
|
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
||||||
|
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
||||||
|
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
||||||
|
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
||||||
|
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
||||||
|
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
||||||
|
*/ |
||||||
|
|
||||||
|
#ifndef ARCH_H |
||||||
|
#define ARCH_H |
||||||
|
|
||||||
|
/* A couple test to catch stupid option combinations */ |
||||||
|
#ifdef FIXED_POINT |
||||||
|
|
||||||
|
#ifdef FLOATING_POINT |
||||||
|
#error You cannot compile as floating point and fixed point at the same time |
||||||
|
#endif |
||||||
|
#ifdef USE_SSE |
||||||
|
#error SSE is only for floating-point |
||||||
|
#endif |
||||||
|
#if defined(ARM4_ASM) + defined(ARM5E_ASM) + defined(BFIN_ASM) > 1 |
||||||
|
#error Make up your mind. What CPU do you have? |
||||||
|
#endif |
||||||
|
#ifdef VORBIS_PSYCHO |
||||||
|
#error Vorbis-psy model currently not implemented in fixed-point |
||||||
|
#endif |
||||||
|
|
||||||
|
#else |
||||||
|
|
||||||
|
#ifndef FLOATING_POINT |
||||||
|
#error You now need to define either FIXED_POINT or FLOATING_POINT |
||||||
|
#endif |
||||||
|
#if defined(ARM4_ASM) || defined(ARM5E_ASM) || defined(BFIN_ASM) |
||||||
|
#error I suppose you can have a [ARM4/ARM5E/Blackfin] that has float instructions? |
||||||
|
#endif |
||||||
|
#ifdef FIXED_DEBUG |
||||||
|
#error "Don't you think enabling fixed-point is a good thing to do if you want to debug that?" |
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
#endif |
||||||
|
|
||||||
|
#ifndef OUTSIDE_SPEEX |
||||||
|
#include "speex/speexdsp_types.h" |
||||||
|
#endif |
||||||
|
|
||||||
|
#define ABS(x) ((x) < 0 ? (-(x)) : (x)) /**< Absolute integer value. */ |
||||||
|
#define ABS16(x) ((x) < 0 ? (-(x)) : (x)) /**< Absolute 16-bit value. */ |
||||||
|
#define MIN16(a,b) ((a) < (b) ? (a) : (b)) /**< Maximum 16-bit value. */ |
||||||
|
#define MAX16(a,b) ((a) > (b) ? (a) : (b)) /**< Maximum 16-bit value. */ |
||||||
|
#define ABS32(x) ((x) < 0 ? (-(x)) : (x)) /**< Absolute 32-bit value. */ |
||||||
|
#define MIN32(a,b) ((a) < (b) ? (a) : (b)) /**< Maximum 32-bit value. */ |
||||||
|
#define MAX32(a,b) ((a) > (b) ? (a) : (b)) /**< Maximum 32-bit value. */ |
||||||
|
|
||||||
|
#ifdef FIXED_POINT |
||||||
|
|
||||||
|
typedef spx_int16_t spx_word16_t; |
||||||
|
typedef spx_int32_t spx_word32_t; |
||||||
|
typedef spx_word32_t spx_mem_t; |
||||||
|
typedef spx_word16_t spx_coef_t; |
||||||
|
typedef spx_word16_t spx_lsp_t; |
||||||
|
typedef spx_word32_t spx_sig_t; |
||||||
|
|
||||||
|
#define Q15ONE 32767 |
||||||
|
|
||||||
|
#define LPC_SCALING 8192 |
||||||
|
#define SIG_SCALING 16384 |
||||||
|
#define LSP_SCALING 8192. |
||||||
|
#define GAMMA_SCALING 32768. |
||||||
|
#define GAIN_SCALING 64 |
||||||
|
#define GAIN_SCALING_1 0.015625 |
||||||
|
|
||||||
|
#define LPC_SHIFT 13 |
||||||
|
#define LSP_SHIFT 13 |
||||||
|
#define SIG_SHIFT 14 |
||||||
|
#define GAIN_SHIFT 6 |
||||||
|
|
||||||
|
#define WORD2INT(x) ((x) < -32767 ? -32768 : ((x) > 32766 ? 32767 : (x))) |
||||||
|
|
||||||
|
#define VERY_SMALL 0 |
||||||
|
#define VERY_LARGE32 ((spx_word32_t)2147483647) |
||||||
|
#define VERY_LARGE16 ((spx_word16_t)32767) |
||||||
|
#define Q15_ONE ((spx_word16_t)32767) |
||||||
|
|
||||||
|
|
||||||
|
#ifdef FIXED_DEBUG |
||||||
|
#include "fixed_debug.h" |
||||||
|
#else |
||||||
|
|
||||||
|
#include "fixed_generic.h" |
||||||
|
|
||||||
|
#ifdef ARM5E_ASM |
||||||
|
#include "fixed_arm5e.h" |
||||||
|
#elif defined(ARM4_ASM) |
||||||
|
#include "fixed_arm4.h" |
||||||
|
#elif defined(BFIN_ASM) |
||||||
|
#include "fixed_bfin.h" |
||||||
|
#endif |
||||||
|
|
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
#else |
||||||
|
|
||||||
|
typedef float spx_mem_t; |
||||||
|
typedef float spx_coef_t; |
||||||
|
typedef float spx_lsp_t; |
||||||
|
typedef float spx_sig_t; |
||||||
|
typedef float spx_word16_t; |
||||||
|
typedef float spx_word32_t; |
||||||
|
|
||||||
|
#define Q15ONE 1.0f |
||||||
|
#define LPC_SCALING 1.f |
||||||
|
#define SIG_SCALING 1.f |
||||||
|
#define LSP_SCALING 1.f |
||||||
|
#define GAMMA_SCALING 1.f |
||||||
|
#define GAIN_SCALING 1.f |
||||||
|
#define GAIN_SCALING_1 1.f |
||||||
|
|
||||||
|
|
||||||
|
#define VERY_SMALL 1e-15f |
||||||
|
#define VERY_LARGE32 1e15f |
||||||
|
#define VERY_LARGE16 1e15f |
||||||
|
#define Q15_ONE ((spx_word16_t)1.f) |
||||||
|
|
||||||
|
#define QCONST16(x,bits) (x) |
||||||
|
#define QCONST32(x,bits) (x) |
||||||
|
|
||||||
|
#define NEG16(x) (-(x)) |
||||||
|
#define NEG32(x) (-(x)) |
||||||
|
#define EXTRACT16(x) (x) |
||||||
|
#define EXTEND32(x) (x) |
||||||
|
#define SHR16(a,shift) (a) |
||||||
|
#define SHL16(a,shift) (a) |
||||||
|
#define SHR32(a,shift) (a) |
||||||
|
#define SHL32(a,shift) (a) |
||||||
|
#define PSHR16(a,shift) (a) |
||||||
|
#define PSHR32(a,shift) (a) |
||||||
|
#define VSHR32(a,shift) (a) |
||||||
|
#define SATURATE16(x,a) (x) |
||||||
|
#define SATURATE32(x,a) (x) |
||||||
|
#define SATURATE32PSHR(x,shift,a) (x) |
||||||
|
|
||||||
|
#define PSHR(a,shift) (a) |
||||||
|
#define SHR(a,shift) (a) |
||||||
|
#define SHL(a,shift) (a) |
||||||
|
#define SATURATE(x,a) (x) |
||||||
|
|
||||||
|
#define ADD16(a,b) ((a)+(b)) |
||||||
|
#define SUB16(a,b) ((a)-(b)) |
||||||
|
#define ADD32(a,b) ((a)+(b)) |
||||||
|
#define SUB32(a,b) ((a)-(b)) |
||||||
|
#define MULT16_16_16(a,b) ((a)*(b)) |
||||||
|
#define MULT16_32_32(a,b) ((a)*(b)) |
||||||
|
#define MULT16_16(a,b) ((spx_word32_t)(a)*(spx_word32_t)(b)) |
||||||
|
#define MAC16_16(c,a,b) ((c)+(spx_word32_t)(a)*(spx_word32_t)(b)) |
||||||
|
|
||||||
|
#define MULT16_32_Q15(a,b) ((a)*(b)) |
||||||
|
#define MULT16_32_P15(a,b) ((a)*(b)) |
||||||
|
|
||||||
|
#define MAC16_32_Q15(c,a,b) ((c)+(a)*(b)) |
||||||
|
|
||||||
|
#define MAC16_16_Q11(c,a,b) ((c)+(a)*(b)) |
||||||
|
#define MAC16_16_Q13(c,a,b) ((c)+(a)*(b)) |
||||||
|
#define MAC16_16_P13(c,a,b) ((c)+(a)*(b)) |
||||||
|
#define MULT16_16_Q11_32(a,b) ((a)*(b)) |
||||||
|
#define MULT16_16_Q13(a,b) ((a)*(b)) |
||||||
|
#define MULT16_16_Q14(a,b) ((a)*(b)) |
||||||
|
#define MULT16_16_Q15(a,b) ((a)*(b)) |
||||||
|
#define MULT16_16_P15(a,b) ((a)*(b)) |
||||||
|
#define MULT16_16_P13(a,b) ((a)*(b)) |
||||||
|
#define MULT16_16_P14(a,b) ((a)*(b)) |
||||||
|
|
||||||
|
#define DIV32_16(a,b) (((spx_word32_t)(a))/(spx_word16_t)(b)) |
||||||
|
#define PDIV32_16(a,b) (((spx_word32_t)(a))/(spx_word16_t)(b)) |
||||||
|
#define DIV32(a,b) (((spx_word32_t)(a))/(spx_word32_t)(b)) |
||||||
|
#define PDIV32(a,b) (((spx_word32_t)(a))/(spx_word32_t)(b)) |
||||||
|
|
||||||
|
#define WORD2INT(x) ((x) < -32767.5f ? -32768 : \ |
||||||
|
((x) > 32766.5f ? 32767 : (spx_int16_t)floor(.5 + (x)))) |
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
#if defined(CONFIG_TI_C54X) || defined(CONFIG_TI_C55X) |
||||||
|
|
||||||
|
/* 2 on TI C5x DSP */ |
||||||
|
#define BYTES_PER_CHAR 2 |
||||||
|
#define BITS_PER_CHAR 16 |
||||||
|
#define LOG2_BITS_PER_CHAR 4 |
||||||
|
|
||||||
|
#else |
||||||
|
|
||||||
|
#define BYTES_PER_CHAR 1 |
||||||
|
#define BITS_PER_CHAR 8 |
||||||
|
#define LOG2_BITS_PER_CHAR 3 |
||||||
|
|
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#ifdef FIXED_DEBUG |
||||||
|
extern long long spx_mips; |
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
#endif |
||||||
@ -0,0 +1,448 @@ |
|||||||
|
/* Copyright (C) 2005-2006 Jean-Marc Valin
|
||||||
|
File: fftwrap.c |
||||||
|
|
||||||
|
Wrapper for various FFTs |
||||||
|
|
||||||
|
Redistribution and use in source and binary forms, with or without |
||||||
|
modification, are permitted provided that the following conditions |
||||||
|
are met: |
||||||
|
|
||||||
|
- Redistributions of source code must retain the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer. |
||||||
|
|
||||||
|
- Redistributions in binary form must reproduce the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer in the |
||||||
|
documentation and/or other materials provided with the distribution. |
||||||
|
|
||||||
|
- Neither the name of the Xiph.org Foundation nor the names of its |
||||||
|
contributors may be used to endorse or promote products derived from |
||||||
|
this software without specific prior written permission. |
||||||
|
|
||||||
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
||||||
|
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
||||||
|
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
||||||
|
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR |
||||||
|
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
||||||
|
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
||||||
|
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
||||||
|
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
||||||
|
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
||||||
|
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
||||||
|
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
||||||
|
|
||||||
|
*/ |
||||||
|
|
||||||
|
#ifdef HAVE_CONFIG_H |
||||||
|
#include "config.h" |
||||||
|
#endif |
||||||
|
|
||||||
|
#include "arch.h" |
||||||
|
#include "os_support.h" |
||||||
|
|
||||||
|
#define MAX_FFT_SIZE 2048 |
||||||
|
|
||||||
|
#ifdef FIXED_POINT |
||||||
|
static int maximize_range(spx_word16_t *in, spx_word16_t *out, spx_word16_t bound, int len) |
||||||
|
{ |
||||||
|
int i, shift; |
||||||
|
spx_word16_t max_val = 0; |
||||||
|
for (i=0;i<len;i++) |
||||||
|
{ |
||||||
|
if (in[i]>max_val) |
||||||
|
max_val = in[i]; |
||||||
|
if (-in[i]>max_val) |
||||||
|
max_val = -in[i]; |
||||||
|
} |
||||||
|
shift=0; |
||||||
|
while (max_val <= (bound>>1) && max_val != 0) |
||||||
|
{ |
||||||
|
max_val <<= 1; |
||||||
|
shift++; |
||||||
|
} |
||||||
|
for (i=0;i<len;i++) |
||||||
|
{ |
||||||
|
out[i] = SHL16(in[i], shift); |
||||||
|
} |
||||||
|
return shift; |
||||||
|
} |
||||||
|
|
||||||
|
static void renorm_range(spx_word16_t *in, spx_word16_t *out, int shift, int len) |
||||||
|
{ |
||||||
|
int i; |
||||||
|
for (i=0;i<len;i++) |
||||||
|
{ |
||||||
|
out[i] = PSHR16(in[i], shift); |
||||||
|
} |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
#ifdef USE_SMALLFT |
||||||
|
|
||||||
|
#include "smallft.h" |
||||||
|
#include <math.h> |
||||||
|
|
||||||
|
void *spx_fft_init(int size) |
||||||
|
{ |
||||||
|
struct drft_lookup *table; |
||||||
|
table = speex_alloc(sizeof(struct drft_lookup)); |
||||||
|
spx_drft_init((struct drft_lookup *)table, size); |
||||||
|
return (void*)table; |
||||||
|
} |
||||||
|
|
||||||
|
void spx_fft_destroy(void *table) |
||||||
|
{ |
||||||
|
spx_drft_clear(table); |
||||||
|
speex_free(table); |
||||||
|
} |
||||||
|
|
||||||
|
void spx_fft(void *table, float *in, float *out) |
||||||
|
{ |
||||||
|
if (in==out) |
||||||
|
{ |
||||||
|
int i; |
||||||
|
float scale = 1./((struct drft_lookup *)table)->n; |
||||||
|
speex_warning("FFT should not be done in-place"); |
||||||
|
for (i=0;i<((struct drft_lookup *)table)->n;i++) |
||||||
|
out[i] = scale*in[i]; |
||||||
|
} else { |
||||||
|
int i; |
||||||
|
float scale = 1./((struct drft_lookup *)table)->n; |
||||||
|
for (i=0;i<((struct drft_lookup *)table)->n;i++) |
||||||
|
out[i] = scale*in[i]; |
||||||
|
} |
||||||
|
spx_drft_forward((struct drft_lookup *)table, out); |
||||||
|
} |
||||||
|
|
||||||
|
void spx_ifft(void *table, float *in, float *out) |
||||||
|
{ |
||||||
|
if (in==out) |
||||||
|
{ |
||||||
|
speex_warning("FFT should not be done in-place"); |
||||||
|
} else { |
||||||
|
int i; |
||||||
|
for (i=0;i<((struct drft_lookup *)table)->n;i++) |
||||||
|
out[i] = in[i]; |
||||||
|
} |
||||||
|
spx_drft_backward((struct drft_lookup *)table, out); |
||||||
|
} |
||||||
|
|
||||||
|
#elif defined(USE_INTEL_MKL) |
||||||
|
#include <mkl.h> |
||||||
|
|
||||||
|
struct mkl_config { |
||||||
|
DFTI_DESCRIPTOR_HANDLE desc; |
||||||
|
int N; |
||||||
|
}; |
||||||
|
|
||||||
|
void *spx_fft_init(int size) |
||||||
|
{ |
||||||
|
struct mkl_config *table = (struct mkl_config *) speex_alloc(sizeof(struct mkl_config)); |
||||||
|
table->N = size; |
||||||
|
DftiCreateDescriptor(&table->desc, DFTI_SINGLE, DFTI_REAL, 1, size); |
||||||
|
DftiSetValue(table->desc, DFTI_PACKED_FORMAT, DFTI_PACK_FORMAT); |
||||||
|
DftiSetValue(table->desc, DFTI_PLACEMENT, DFTI_NOT_INPLACE); |
||||||
|
DftiSetValue(table->desc, DFTI_FORWARD_SCALE, 1.0f / size); |
||||||
|
DftiCommitDescriptor(table->desc); |
||||||
|
return table; |
||||||
|
} |
||||||
|
|
||||||
|
void spx_fft_destroy(void *table) |
||||||
|
{ |
||||||
|
struct mkl_config *t = (struct mkl_config *) table; |
||||||
|
DftiFreeDescriptor(t->desc); |
||||||
|
speex_free(table); |
||||||
|
} |
||||||
|
|
||||||
|
void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) |
||||||
|
{ |
||||||
|
struct mkl_config *t = (struct mkl_config *) table; |
||||||
|
DftiComputeForward(t->desc, in, out); |
||||||
|
} |
||||||
|
|
||||||
|
void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out) |
||||||
|
{ |
||||||
|
struct mkl_config *t = (struct mkl_config *) table; |
||||||
|
DftiComputeBackward(t->desc, in, out); |
||||||
|
} |
||||||
|
|
||||||
|
#elif defined(USE_INTEL_IPP) |
||||||
|
|
||||||
|
#include <ipps.h> |
||||||
|
|
||||||
|
struct ipp_fft_config |
||||||
|
{ |
||||||
|
IppsDFTSpec_R_32f *dftSpec; |
||||||
|
Ipp8u *buffer; |
||||||
|
}; |
||||||
|
|
||||||
|
void *spx_fft_init(int size) |
||||||
|
{ |
||||||
|
int bufferSize = 0; |
||||||
|
int hint; |
||||||
|
struct ipp_fft_config *table; |
||||||
|
|
||||||
|
table = (struct ipp_fft_config *)speex_alloc(sizeof(struct ipp_fft_config)); |
||||||
|
|
||||||
|
/* there appears to be no performance difference between ippAlgHintFast and
|
||||||
|
ippAlgHintAccurate when using the with the floating point version |
||||||
|
of the fft. */ |
||||||
|
hint = ippAlgHintAccurate; |
||||||
|
|
||||||
|
ippsDFTInitAlloc_R_32f(&table->dftSpec, size, IPP_FFT_DIV_FWD_BY_N, hint); |
||||||
|
|
||||||
|
ippsDFTGetBufSize_R_32f(table->dftSpec, &bufferSize); |
||||||
|
table->buffer = ippsMalloc_8u(bufferSize); |
||||||
|
|
||||||
|
return table; |
||||||
|
} |
||||||
|
|
||||||
|
void spx_fft_destroy(void *table) |
||||||
|
{ |
||||||
|
struct ipp_fft_config *t = (struct ipp_fft_config *)table; |
||||||
|
ippsFree(t->buffer); |
||||||
|
ippsDFTFree_R_32f(t->dftSpec); |
||||||
|
speex_free(t); |
||||||
|
} |
||||||
|
|
||||||
|
void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) |
||||||
|
{ |
||||||
|
struct ipp_fft_config *t = (struct ipp_fft_config *)table; |
||||||
|
ippsDFTFwd_RToPack_32f(in, out, t->dftSpec, t->buffer); |
||||||
|
} |
||||||
|
|
||||||
|
void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out) |
||||||
|
{ |
||||||
|
struct ipp_fft_config *t = (struct ipp_fft_config *)table; |
||||||
|
ippsDFTInv_PackToR_32f(in, out, t->dftSpec, t->buffer); |
||||||
|
} |
||||||
|
|
||||||
|
#elif defined(USE_GPL_FFTW3) |
||||||
|
|
||||||
|
#include <fftw3.h> |
||||||
|
|
||||||
|
struct fftw_config { |
||||||
|
float *in; |
||||||
|
float *out; |
||||||
|
fftwf_plan fft; |
||||||
|
fftwf_plan ifft; |
||||||
|
int N; |
||||||
|
}; |
||||||
|
|
||||||
|
void *spx_fft_init(int size) |
||||||
|
{ |
||||||
|
struct fftw_config *table = (struct fftw_config *) speex_alloc(sizeof(struct fftw_config)); |
||||||
|
table->in = fftwf_malloc(sizeof(float) * (size+2)); |
||||||
|
table->out = fftwf_malloc(sizeof(float) * (size+2)); |
||||||
|
|
||||||
|
table->fft = fftwf_plan_dft_r2c_1d(size, table->in, (fftwf_complex *) table->out, FFTW_PATIENT); |
||||||
|
table->ifft = fftwf_plan_dft_c2r_1d(size, (fftwf_complex *) table->in, table->out, FFTW_PATIENT); |
||||||
|
|
||||||
|
table->N = size; |
||||||
|
return table; |
||||||
|
} |
||||||
|
|
||||||
|
void spx_fft_destroy(void *table) |
||||||
|
{ |
||||||
|
struct fftw_config *t = (struct fftw_config *) table; |
||||||
|
fftwf_destroy_plan(t->fft); |
||||||
|
fftwf_destroy_plan(t->ifft); |
||||||
|
fftwf_free(t->in); |
||||||
|
fftwf_free(t->out); |
||||||
|
speex_free(table); |
||||||
|
} |
||||||
|
|
||||||
|
|
||||||
|
void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) |
||||||
|
{ |
||||||
|
int i; |
||||||
|
struct fftw_config *t = (struct fftw_config *) table; |
||||||
|
const int N = t->N; |
||||||
|
float *iptr = t->in; |
||||||
|
float *optr = t->out; |
||||||
|
const float m = 1.0 / N; |
||||||
|
for(i=0;i<N;++i) |
||||||
|
iptr[i]=in[i] * m; |
||||||
|
|
||||||
|
fftwf_execute(t->fft); |
||||||
|
|
||||||
|
out[0] = optr[0]; |
||||||
|
for(i=1;i<N;++i) |
||||||
|
out[i] = optr[i+1]; |
||||||
|
} |
||||||
|
|
||||||
|
void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out) |
||||||
|
{ |
||||||
|
int i; |
||||||
|
struct fftw_config *t = (struct fftw_config *) table; |
||||||
|
const int N = t->N; |
||||||
|
float *iptr = t->in; |
||||||
|
float *optr = t->out; |
||||||
|
|
||||||
|
iptr[0] = in[0]; |
||||||
|
iptr[1] = 0.0f; |
||||||
|
for(i=1;i<N;++i) |
||||||
|
iptr[i+1] = in[i]; |
||||||
|
iptr[N+1] = 0.0f; |
||||||
|
|
||||||
|
fftwf_execute(t->ifft); |
||||||
|
|
||||||
|
for(i=0;i<N;++i) |
||||||
|
out[i] = optr[i]; |
||||||
|
} |
||||||
|
|
||||||
|
#elif defined(USE_KISS_FFT) |
||||||
|
|
||||||
|
#include "kiss_fftr.h" |
||||||
|
#include "kiss_fft.h" |
||||||
|
|
||||||
|
struct kiss_config { |
||||||
|
kiss_fftr_cfg forward; |
||||||
|
kiss_fftr_cfg backward; |
||||||
|
int N; |
||||||
|
}; |
||||||
|
|
||||||
|
void *spx_fft_init(int size) |
||||||
|
{ |
||||||
|
struct kiss_config *table; |
||||||
|
table = (struct kiss_config*)speex_alloc(sizeof(struct kiss_config)); |
||||||
|
table->forward = kiss_fftr_alloc(size,0,NULL,NULL); |
||||||
|
table->backward = kiss_fftr_alloc(size,1,NULL,NULL); |
||||||
|
table->N = size; |
||||||
|
return table; |
||||||
|
} |
||||||
|
|
||||||
|
void spx_fft_destroy(void *table) |
||||||
|
{ |
||||||
|
struct kiss_config *t = (struct kiss_config *)table; |
||||||
|
kiss_fftr_free(t->forward); |
||||||
|
kiss_fftr_free(t->backward); |
||||||
|
speex_free(table); |
||||||
|
} |
||||||
|
|
||||||
|
#ifdef FIXED_POINT |
||||||
|
|
||||||
|
void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) |
||||||
|
{ |
||||||
|
int shift; |
||||||
|
struct kiss_config *t = (struct kiss_config *)table; |
||||||
|
shift = maximize_range(in, in, 32000, t->N); |
||||||
|
kiss_fftr2(t->forward, in, out); |
||||||
|
renorm_range(in, in, shift, t->N); |
||||||
|
renorm_range(out, out, shift, t->N); |
||||||
|
} |
||||||
|
|
||||||
|
#else |
||||||
|
|
||||||
|
void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out) |
||||||
|
{ |
||||||
|
int i; |
||||||
|
float scale; |
||||||
|
struct kiss_config *t = (struct kiss_config *)table; |
||||||
|
scale = 1./t->N; |
||||||
|
kiss_fftr2(t->forward, in, out); |
||||||
|
for (i=0;i<t->N;i++) |
||||||
|
out[i] *= scale; |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out) |
||||||
|
{ |
||||||
|
struct kiss_config *t = (struct kiss_config *)table; |
||||||
|
kiss_fftri2(t->backward, in, out); |
||||||
|
} |
||||||
|
|
||||||
|
|
||||||
|
#else |
||||||
|
|
||||||
|
#error No other FFT implemented |
||||||
|
|
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
#ifdef FIXED_POINT |
||||||
|
/*#include "smallft.h"*/ |
||||||
|
|
||||||
|
|
||||||
|
void spx_fft_float(void *table, float *in, float *out) |
||||||
|
{ |
||||||
|
int i; |
||||||
|
#ifdef USE_SMALLFT |
||||||
|
int N = ((struct drft_lookup *)table)->n; |
||||||
|
#elif defined(USE_KISS_FFT) |
||||||
|
int N = ((struct kiss_config *)table)->N; |
||||||
|
#else |
||||||
|
#endif |
||||||
|
#ifdef VAR_ARRAYS |
||||||
|
spx_word16_t _in[N]; |
||||||
|
spx_word16_t _out[N]; |
||||||
|
#else |
||||||
|
spx_word16_t _in[MAX_FFT_SIZE]; |
||||||
|
spx_word16_t _out[MAX_FFT_SIZE]; |
||||||
|
#endif |
||||||
|
for (i=0;i<N;i++) |
||||||
|
_in[i] = (int)floor(.5+in[i]); |
||||||
|
spx_fft(table, _in, _out); |
||||||
|
for (i=0;i<N;i++) |
||||||
|
out[i] = _out[i]; |
||||||
|
#if 0 |
||||||
|
if (!fixed_point) |
||||||
|
{ |
||||||
|
float scale; |
||||||
|
struct drft_lookup t; |
||||||
|
spx_drft_init(&t, ((struct kiss_config *)table)->N); |
||||||
|
scale = 1./((struct kiss_config *)table)->N; |
||||||
|
for (i=0;i<((struct kiss_config *)table)->N;i++) |
||||||
|
out[i] = scale*in[i]; |
||||||
|
spx_drft_forward(&t, out); |
||||||
|
spx_drft_clear(&t); |
||||||
|
} |
||||||
|
#endif |
||||||
|
} |
||||||
|
|
||||||
|
void spx_ifft_float(void *table, float *in, float *out) |
||||||
|
{ |
||||||
|
int i; |
||||||
|
#ifdef USE_SMALLFT |
||||||
|
int N = ((struct drft_lookup *)table)->n; |
||||||
|
#elif defined(USE_KISS_FFT) |
||||||
|
int N = ((struct kiss_config *)table)->N; |
||||||
|
#else |
||||||
|
#endif |
||||||
|
#ifdef VAR_ARRAYS |
||||||
|
spx_word16_t _in[N]; |
||||||
|
spx_word16_t _out[N]; |
||||||
|
#else |
||||||
|
spx_word16_t _in[MAX_FFT_SIZE]; |
||||||
|
spx_word16_t _out[MAX_FFT_SIZE]; |
||||||
|
#endif |
||||||
|
for (i=0;i<N;i++) |
||||||
|
_in[i] = (int)floor(.5+in[i]); |
||||||
|
spx_ifft(table, _in, _out); |
||||||
|
for (i=0;i<N;i++) |
||||||
|
out[i] = _out[i]; |
||||||
|
#if 0 |
||||||
|
if (!fixed_point) |
||||||
|
{ |
||||||
|
int i; |
||||||
|
struct drft_lookup t; |
||||||
|
spx_drft_init(&t, ((struct kiss_config *)table)->N); |
||||||
|
for (i=0;i<((struct kiss_config *)table)->N;i++) |
||||||
|
out[i] = in[i]; |
||||||
|
spx_drft_backward(&t, out); |
||||||
|
spx_drft_clear(&t); |
||||||
|
} |
||||||
|
#endif |
||||||
|
} |
||||||
|
|
||||||
|
#else |
||||||
|
|
||||||
|
void spx_fft_float(void *table, float *in, float *out) |
||||||
|
{ |
||||||
|
spx_fft(table, in, out); |
||||||
|
} |
||||||
|
void spx_ifft_float(void *table, float *in, float *out) |
||||||
|
{ |
||||||
|
spx_ifft(table, in, out); |
||||||
|
} |
||||||
|
|
||||||
|
#endif |
||||||
@ -0,0 +1,58 @@ |
|||||||
|
/* Copyright (C) 2005 Jean-Marc Valin
|
||||||
|
File: fftwrap.h |
||||||
|
|
||||||
|
Wrapper for various FFTs |
||||||
|
|
||||||
|
Redistribution and use in source and binary forms, with or without |
||||||
|
modification, are permitted provided that the following conditions |
||||||
|
are met: |
||||||
|
|
||||||
|
- Redistributions of source code must retain the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer. |
||||||
|
|
||||||
|
- Redistributions in binary form must reproduce the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer in the |
||||||
|
documentation and/or other materials provided with the distribution. |
||||||
|
|
||||||
|
- Neither the name of the Xiph.org Foundation nor the names of its |
||||||
|
contributors may be used to endorse or promote products derived from |
||||||
|
this software without specific prior written permission. |
||||||
|
|
||||||
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
||||||
|
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
||||||
|
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
||||||
|
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR |
||||||
|
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
||||||
|
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
||||||
|
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
||||||
|
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
||||||
|
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
||||||
|
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
||||||
|
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
||||||
|
|
||||||
|
*/ |
||||||
|
|
||||||
|
#ifndef FFTWRAP_H |
||||||
|
#define FFTWRAP_H |
||||||
|
|
||||||
|
#include "arch.h" |
||||||
|
|
||||||
|
/** Compute tables for an FFT */ |
||||||
|
void *spx_fft_init(int size); |
||||||
|
|
||||||
|
/** Destroy tables for an FFT */ |
||||||
|
void spx_fft_destroy(void *table); |
||||||
|
|
||||||
|
/** Forward (real to half-complex) transform */ |
||||||
|
void spx_fft(void *table, spx_word16_t *in, spx_word16_t *out); |
||||||
|
|
||||||
|
/** Backward (half-complex to real) transform */ |
||||||
|
void spx_ifft(void *table, spx_word16_t *in, spx_word16_t *out); |
||||||
|
|
||||||
|
/** Forward (real to half-complex) transform of float data */ |
||||||
|
void spx_fft_float(void *table, float *in, float *out); |
||||||
|
|
||||||
|
/** Backward (half-complex to real) transform of float data */ |
||||||
|
void spx_ifft_float(void *table, float *in, float *out); |
||||||
|
|
||||||
|
#endif |
||||||
@ -0,0 +1,170 @@ |
|||||||
|
/* Copyright (C) Jean-Marc Valin */ |
||||||
|
/**
|
||||||
|
@file speex_echo.h |
||||||
|
@brief Echo cancellation |
||||||
|
*/ |
||||||
|
/*
|
||||||
|
Redistribution and use in source and binary forms, with or without |
||||||
|
modification, are permitted provided that the following conditions are |
||||||
|
met: |
||||||
|
|
||||||
|
1. Redistributions of source code must retain the above copyright notice, |
||||||
|
this list of conditions and the following disclaimer. |
||||||
|
|
||||||
|
2. Redistributions in binary form must reproduce the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer in the |
||||||
|
documentation and/or other materials provided with the distribution. |
||||||
|
|
||||||
|
3. The name of the author may not be used to endorse or promote products |
||||||
|
derived from this software without specific prior written permission. |
||||||
|
|
||||||
|
THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR |
||||||
|
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES |
||||||
|
OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE |
||||||
|
DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, |
||||||
|
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES |
||||||
|
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR |
||||||
|
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
||||||
|
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, |
||||||
|
STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN |
||||||
|
ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
||||||
|
POSSIBILITY OF SUCH DAMAGE. |
||||||
|
*/ |
||||||
|
|
||||||
|
#ifndef SPEEX_ECHO_H |
||||||
|
#define SPEEX_ECHO_H |
||||||
|
/** @defgroup SpeexEchoState SpeexEchoState: Acoustic echo canceller
|
||||||
|
* This is the acoustic echo canceller module. |
||||||
|
* @{ |
||||||
|
*/ |
||||||
|
#include "speexdsp_types.h" |
||||||
|
|
||||||
|
#ifdef __cplusplus |
||||||
|
extern "C" { |
||||||
|
#endif |
||||||
|
|
||||||
|
/** Obtain frame size used by the AEC */ |
||||||
|
#define SPEEX_ECHO_GET_FRAME_SIZE 3 |
||||||
|
|
||||||
|
/** Set sampling rate */ |
||||||
|
#define SPEEX_ECHO_SET_SAMPLING_RATE 24 |
||||||
|
/** Get sampling rate */ |
||||||
|
#define SPEEX_ECHO_GET_SAMPLING_RATE 25 |
||||||
|
|
||||||
|
/* Can't set window sizes */ |
||||||
|
/** Get size of impulse response (int32) */ |
||||||
|
#define SPEEX_ECHO_GET_IMPULSE_RESPONSE_SIZE 27 |
||||||
|
|
||||||
|
/* Can't set window content */ |
||||||
|
/** Get impulse response (int32[]) */ |
||||||
|
#define SPEEX_ECHO_GET_IMPULSE_RESPONSE 29 |
||||||
|
|
||||||
|
/** Internal echo canceller state. Should never be accessed directly. */ |
||||||
|
struct SpeexEchoState_; |
||||||
|
|
||||||
|
/** @class SpeexEchoState
|
||||||
|
* This holds the state of the echo canceller. You need one per channel. |
||||||
|
*/ |
||||||
|
|
||||||
|
/** Internal echo canceller state. Should never be accessed directly. */ |
||||||
|
typedef struct SpeexEchoState_ SpeexEchoState; |
||||||
|
|
||||||
|
/** Creates a new echo canceller state
|
||||||
|
* @param frame_size Number of samples to process at one time (should correspond to 10-20 ms) |
||||||
|
* @param filter_length Number of samples of echo to cancel (should generally correspond to 100-500 ms) |
||||||
|
* @return Newly-created echo canceller state |
||||||
|
*/ |
||||||
|
SpeexEchoState *speex_echo_state_init(int frame_size, int filter_length); |
||||||
|
|
||||||
|
/** Creates a new multi-channel echo canceller state
|
||||||
|
* @param frame_size Number of samples to process at one time (should correspond to 10-20 ms) |
||||||
|
* @param filter_length Number of samples of echo to cancel (should generally correspond to 100-500 ms) |
||||||
|
* @param nb_mic Number of microphone channels |
||||||
|
* @param nb_speakers Number of speaker channels |
||||||
|
* @return Newly-created echo canceller state |
||||||
|
*/ |
||||||
|
SpeexEchoState *speex_echo_state_init_mc(int frame_size, int filter_length, int nb_mic, int nb_speakers); |
||||||
|
|
||||||
|
/** Destroys an echo canceller state
|
||||||
|
* @param st Echo canceller state |
||||||
|
*/ |
||||||
|
void speex_echo_state_destroy(SpeexEchoState *st); |
||||||
|
|
||||||
|
/** Performs echo cancellation a frame, based on the audio sent to the speaker (no delay is added
|
||||||
|
* to playback in this form) |
||||||
|
* |
||||||
|
* @param st Echo canceller state |
||||||
|
* @param rec Signal from the microphone (near end + far end echo) |
||||||
|
* @param play Signal played to the speaker (received from far end) |
||||||
|
* @param out Returns near-end signal with echo removed |
||||||
|
*/ |
||||||
|
void speex_echo_cancellation(SpeexEchoState *st, const spx_int16_t *rec, const spx_int16_t *play, spx_int16_t *out); |
||||||
|
|
||||||
|
/** Performs echo cancellation a frame (deprecated) */ |
||||||
|
void speex_echo_cancel(SpeexEchoState *st, const spx_int16_t *rec, const spx_int16_t *play, spx_int16_t *out, spx_int32_t *Yout); |
||||||
|
|
||||||
|
/** Perform echo cancellation using internal playback buffer, which is delayed by two frames
|
||||||
|
* to account for the delay introduced by most soundcards (but it could be off!) |
||||||
|
* @param st Echo canceller state |
||||||
|
* @param rec Signal from the microphone (near end + far end echo) |
||||||
|
* @param out Returns near-end signal with echo removed |
||||||
|
*/ |
||||||
|
void speex_echo_capture(SpeexEchoState *st, const spx_int16_t *rec, spx_int16_t *out); |
||||||
|
|
||||||
|
/** Let the echo canceller know that a frame was just queued to the soundcard
|
||||||
|
* @param st Echo canceller state |
||||||
|
* @param play Signal played to the speaker (received from far end) |
||||||
|
*/ |
||||||
|
void speex_echo_playback(SpeexEchoState *st, const spx_int16_t *play); |
||||||
|
|
||||||
|
/** Reset the echo canceller to its original state
|
||||||
|
* @param st Echo canceller state |
||||||
|
*/ |
||||||
|
void speex_echo_state_reset(SpeexEchoState *st); |
||||||
|
|
||||||
|
/** Used like the ioctl function to control the echo canceller parameters
|
||||||
|
* |
||||||
|
* @param st Echo canceller state |
||||||
|
* @param request ioctl-type request (one of the SPEEX_ECHO_* macros) |
||||||
|
* @param ptr Data exchanged to-from function |
||||||
|
* @return 0 if no error, -1 if request in unknown |
||||||
|
*/ |
||||||
|
int speex_echo_ctl(SpeexEchoState *st, int request, void *ptr); |
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
struct SpeexDecorrState_; |
||||||
|
|
||||||
|
typedef struct SpeexDecorrState_ SpeexDecorrState; |
||||||
|
|
||||||
|
|
||||||
|
/** Create a state for the channel decorrelation algorithm
|
||||||
|
This is useful for multi-channel echo cancellation only |
||||||
|
* @param rate Sampling rate |
||||||
|
* @param channels Number of channels (it's a bit pointless if you don't have at least 2) |
||||||
|
* @param frame_size Size of the frame to process at ones (counting samples *per* channel) |
||||||
|
*/ |
||||||
|
SpeexDecorrState *speex_decorrelate_new(int rate, int channels, int frame_size); |
||||||
|
|
||||||
|
/** Remove correlation between the channels by modifying the phase and possibly
|
||||||
|
adding noise in a way that is not (or little) perceptible. |
||||||
|
* @param st Decorrelator state |
||||||
|
* @param in Input audio in interleaved format |
||||||
|
* @param out Result of the decorrelation (out *may* alias in) |
||||||
|
* @param strength How much alteration of the audio to apply from 0 to 100. |
||||||
|
*/ |
||||||
|
void speex_decorrelate(SpeexDecorrState *st, const spx_int16_t *in, spx_int16_t *out, int strength); |
||||||
|
|
||||||
|
/** Destroy a Decorrelation state
|
||||||
|
* @param st State to destroy |
||||||
|
*/ |
||||||
|
void speex_decorrelate_destroy(SpeexDecorrState *st); |
||||||
|
|
||||||
|
|
||||||
|
#ifdef __cplusplus |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
/** @}*/ |
||||||
|
#endif |
||||||
@ -0,0 +1,14 @@ |
|||||||
|
/* speexdsp_config_types.h — статически сгенерировано для POSIX (Linux/FreeBSD/Android).
|
||||||
|
* В оригинале генерируется configure из speexdsp_config_types.h.in. На Windows ветка |
||||||
|
* _WIN32 в speexdsp_types.h использует нативные типы и этот файл не инклудится. */ |
||||||
|
#ifndef __SPEEX_TYPES_H__ |
||||||
|
#define __SPEEX_TYPES_H__ |
||||||
|
|
||||||
|
#include <stdint.h> |
||||||
|
|
||||||
|
typedef int16_t spx_int16_t; |
||||||
|
typedef uint16_t spx_uint16_t; |
||||||
|
typedef int32_t spx_int32_t; |
||||||
|
typedef uint32_t spx_uint32_t; |
||||||
|
|
||||||
|
#endif |
||||||
@ -0,0 +1,126 @@ |
|||||||
|
/* speexdsp_types.h taken from libogg */ |
||||||
|
/********************************************************************
|
||||||
|
* * |
||||||
|
* THIS FILE IS PART OF THE OggVorbis SOFTWARE CODEC SOURCE CODE. * |
||||||
|
* USE, DISTRIBUTION AND REPRODUCTION OF THIS LIBRARY SOURCE IS * |
||||||
|
* GOVERNED BY A BSD-STYLE SOURCE LICENSE INCLUDED WITH THIS SOURCE * |
||||||
|
* IN 'COPYING'. PLEASE READ THESE TERMS BEFORE DISTRIBUTING. * |
||||||
|
* * |
||||||
|
* THE OggVorbis SOURCE CODE IS (C) COPYRIGHT 1994-2002 * |
||||||
|
* by the Xiph.Org Foundation http://www.xiph.org/ *
|
||||||
|
* * |
||||||
|
******************************************************************** |
||||||
|
|
||||||
|
function: #ifdef jail to whip a few platforms into the UNIX ideal. |
||||||
|
last mod: $Id: os_types.h 7524 2004-08-11 04:20:36Z conrad $ |
||||||
|
|
||||||
|
********************************************************************/ |
||||||
|
/**
|
||||||
|
@file speexdsp_types.h |
||||||
|
@brief Speex types |
||||||
|
*/ |
||||||
|
#ifndef _SPEEX_TYPES_H |
||||||
|
#define _SPEEX_TYPES_H |
||||||
|
|
||||||
|
#if defined(_WIN32) |
||||||
|
|
||||||
|
# if defined(__CYGWIN__) |
||||||
|
# include <_G_config.h> |
||||||
|
typedef _G_int32_t spx_int32_t; |
||||||
|
typedef _G_uint32_t spx_uint32_t; |
||||||
|
typedef _G_int16_t spx_int16_t; |
||||||
|
typedef _G_uint16_t spx_uint16_t; |
||||||
|
# elif defined(__MINGW32__) |
||||||
|
typedef short spx_int16_t; |
||||||
|
typedef unsigned short spx_uint16_t; |
||||||
|
typedef int spx_int32_t; |
||||||
|
typedef unsigned int spx_uint32_t; |
||||||
|
# elif defined(__MWERKS__) |
||||||
|
typedef int spx_int32_t; |
||||||
|
typedef unsigned int spx_uint32_t; |
||||||
|
typedef short spx_int16_t; |
||||||
|
typedef unsigned short spx_uint16_t; |
||||||
|
# else |
||||||
|
/* MSVC/Borland */ |
||||||
|
typedef __int32 spx_int32_t; |
||||||
|
typedef unsigned __int32 spx_uint32_t; |
||||||
|
typedef __int16 spx_int16_t; |
||||||
|
typedef unsigned __int16 spx_uint16_t; |
||||||
|
# endif |
||||||
|
|
||||||
|
#elif defined(__MACOS__) |
||||||
|
|
||||||
|
# include <sys/types.h> |
||||||
|
typedef SInt16 spx_int16_t; |
||||||
|
typedef UInt16 spx_uint16_t; |
||||||
|
typedef SInt32 spx_int32_t; |
||||||
|
typedef UInt32 spx_uint32_t; |
||||||
|
|
||||||
|
#elif (defined(__APPLE__) && defined(__MACH__)) /* MacOS X Framework build */ |
||||||
|
|
||||||
|
# include <sys/types.h> |
||||||
|
typedef int16_t spx_int16_t; |
||||||
|
typedef u_int16_t spx_uint16_t; |
||||||
|
typedef int32_t spx_int32_t; |
||||||
|
typedef u_int32_t spx_uint32_t; |
||||||
|
|
||||||
|
#elif defined(__BEOS__) |
||||||
|
|
||||||
|
/* Be */ |
||||||
|
# include <inttypes.h> |
||||||
|
typedef int16_t spx_int16_t; |
||||||
|
typedef u_int16_t spx_uint16_t; |
||||||
|
typedef int32_t spx_int32_t; |
||||||
|
typedef u_int32_t spx_uint32_t; |
||||||
|
|
||||||
|
#elif defined (__EMX__) |
||||||
|
|
||||||
|
/* OS/2 GCC */ |
||||||
|
typedef short spx_int16_t; |
||||||
|
typedef unsigned short spx_uint16_t; |
||||||
|
typedef int spx_int32_t; |
||||||
|
typedef unsigned int spx_uint32_t; |
||||||
|
|
||||||
|
#elif defined (DJGPP) |
||||||
|
|
||||||
|
/* DJGPP */ |
||||||
|
typedef short spx_int16_t; |
||||||
|
typedef int spx_int32_t; |
||||||
|
typedef unsigned int spx_uint32_t; |
||||||
|
|
||||||
|
#elif defined(R5900) |
||||||
|
|
||||||
|
/* PS2 EE */ |
||||||
|
typedef int spx_int32_t; |
||||||
|
typedef unsigned spx_uint32_t; |
||||||
|
typedef short spx_int16_t; |
||||||
|
|
||||||
|
#elif defined(__SYMBIAN32__) |
||||||
|
|
||||||
|
/* Symbian GCC */ |
||||||
|
typedef signed short spx_int16_t; |
||||||
|
typedef unsigned short spx_uint16_t; |
||||||
|
typedef signed int spx_int32_t; |
||||||
|
typedef unsigned int spx_uint32_t; |
||||||
|
|
||||||
|
#elif defined(CONFIG_TI_C54X) || defined (CONFIG_TI_C55X) |
||||||
|
|
||||||
|
typedef short spx_int16_t; |
||||||
|
typedef unsigned short spx_uint16_t; |
||||||
|
typedef long spx_int32_t; |
||||||
|
typedef unsigned long spx_uint32_t; |
||||||
|
|
||||||
|
#elif defined(CONFIG_TI_C6X) |
||||||
|
|
||||||
|
typedef short spx_int16_t; |
||||||
|
typedef unsigned short spx_uint16_t; |
||||||
|
typedef int spx_int32_t; |
||||||
|
typedef unsigned int spx_uint32_t; |
||||||
|
|
||||||
|
#else |
||||||
|
|
||||||
|
#include "speexdsp_config_types.h" |
||||||
|
|
||||||
|
#endif |
||||||
|
|
||||||
|
#endif /* _SPEEX_TYPES_H */ |
||||||
@ -0,0 +1,523 @@ |
|||||||
|
/*
|
||||||
|
Copyright (c) 2003-2004, Mark Borgerding |
||||||
|
Copyright (c) 2005-2007, Jean-Marc Valin |
||||||
|
|
||||||
|
All rights reserved. |
||||||
|
|
||||||
|
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: |
||||||
|
|
||||||
|
* Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. |
||||||
|
* Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. |
||||||
|
* Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission. |
||||||
|
|
||||||
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
||||||
|
*/ |
||||||
|
|
||||||
|
|
||||||
|
#ifdef HAVE_CONFIG_H |
||||||
|
#include "config.h" |
||||||
|
#endif |
||||||
|
|
||||||
|
#include "_kiss_fft_guts.h" |
||||||
|
#include "arch.h" |
||||||
|
#include "os_support.h" |
||||||
|
|
||||||
|
/* The guts header contains all the multiplication and addition macros that are defined for
|
||||||
|
fixed or floating point complex numbers. It also declares the kf_ internal functions. |
||||||
|
*/ |
||||||
|
|
||||||
|
static void kf_bfly2( |
||||||
|
kiss_fft_cpx * Fout, |
||||||
|
const size_t fstride, |
||||||
|
const kiss_fft_cfg st, |
||||||
|
int m, |
||||||
|
int N, |
||||||
|
int mm |
||||||
|
) |
||||||
|
{ |
||||||
|
kiss_fft_cpx * Fout2; |
||||||
|
kiss_fft_cpx * tw1; |
||||||
|
kiss_fft_cpx t; |
||||||
|
if (!st->inverse) { |
||||||
|
int i,j; |
||||||
|
kiss_fft_cpx * Fout_beg = Fout; |
||||||
|
for (i=0;i<N;i++) |
||||||
|
{ |
||||||
|
Fout = Fout_beg + i*mm; |
||||||
|
Fout2 = Fout + m; |
||||||
|
tw1 = st->twiddles; |
||||||
|
for(j=0;j<m;j++) |
||||||
|
{ |
||||||
|
/* Almost the same as the code path below, except that we divide the input by two
|
||||||
|
(while keeping the best accuracy possible) */ |
||||||
|
spx_word32_t tr, ti; |
||||||
|
tr = SHR32(SUB32(MULT16_16(Fout2->r , tw1->r),MULT16_16(Fout2->i , tw1->i)), 1); |
||||||
|
ti = SHR32(ADD32(MULT16_16(Fout2->i , tw1->r),MULT16_16(Fout2->r , tw1->i)), 1); |
||||||
|
tw1 += fstride; |
||||||
|
Fout2->r = PSHR32(SUB32(SHL32(EXTEND32(Fout->r), 14), tr), 15); |
||||||
|
Fout2->i = PSHR32(SUB32(SHL32(EXTEND32(Fout->i), 14), ti), 15); |
||||||
|
Fout->r = PSHR32(ADD32(SHL32(EXTEND32(Fout->r), 14), tr), 15); |
||||||
|
Fout->i = PSHR32(ADD32(SHL32(EXTEND32(Fout->i), 14), ti), 15); |
||||||
|
++Fout2; |
||||||
|
++Fout; |
||||||
|
} |
||||||
|
} |
||||||
|
} else { |
||||||
|
int i,j; |
||||||
|
kiss_fft_cpx * Fout_beg = Fout; |
||||||
|
for (i=0;i<N;i++) |
||||||
|
{ |
||||||
|
Fout = Fout_beg + i*mm; |
||||||
|
Fout2 = Fout + m; |
||||||
|
tw1 = st->twiddles; |
||||||
|
for(j=0;j<m;j++) |
||||||
|
{ |
||||||
|
C_MUL (t, *Fout2 , *tw1); |
||||||
|
tw1 += fstride; |
||||||
|
C_SUB( *Fout2 , *Fout , t ); |
||||||
|
C_ADDTO( *Fout , t ); |
||||||
|
++Fout2; |
||||||
|
++Fout; |
||||||
|
} |
||||||
|
} |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
static void kf_bfly4( |
||||||
|
kiss_fft_cpx * Fout, |
||||||
|
const size_t fstride, |
||||||
|
const kiss_fft_cfg st, |
||||||
|
int m, |
||||||
|
int N, |
||||||
|
int mm |
||||||
|
) |
||||||
|
{ |
||||||
|
kiss_fft_cpx *tw1,*tw2,*tw3; |
||||||
|
kiss_fft_cpx scratch[6]; |
||||||
|
const size_t m2=2*m; |
||||||
|
const size_t m3=3*m; |
||||||
|
int i, j; |
||||||
|
|
||||||
|
if (st->inverse) |
||||||
|
{ |
||||||
|
kiss_fft_cpx * Fout_beg = Fout; |
||||||
|
for (i=0;i<N;i++) |
||||||
|
{ |
||||||
|
Fout = Fout_beg + i*mm; |
||||||
|
tw3 = tw2 = tw1 = st->twiddles; |
||||||
|
for (j=0;j<m;j++) |
||||||
|
{ |
||||||
|
C_MUL(scratch[0],Fout[m] , *tw1 ); |
||||||
|
C_MUL(scratch[1],Fout[m2] , *tw2 ); |
||||||
|
C_MUL(scratch[2],Fout[m3] , *tw3 ); |
||||||
|
|
||||||
|
C_SUB( scratch[5] , *Fout, scratch[1] ); |
||||||
|
C_ADDTO(*Fout, scratch[1]); |
||||||
|
C_ADD( scratch[3] , scratch[0] , scratch[2] ); |
||||||
|
C_SUB( scratch[4] , scratch[0] , scratch[2] ); |
||||||
|
C_SUB( Fout[m2], *Fout, scratch[3] ); |
||||||
|
tw1 += fstride; |
||||||
|
tw2 += fstride*2; |
||||||
|
tw3 += fstride*3; |
||||||
|
C_ADDTO( *Fout , scratch[3] ); |
||||||
|
|
||||||
|
Fout[m].r = scratch[5].r - scratch[4].i; |
||||||
|
Fout[m].i = scratch[5].i + scratch[4].r; |
||||||
|
Fout[m3].r = scratch[5].r + scratch[4].i; |
||||||
|
Fout[m3].i = scratch[5].i - scratch[4].r; |
||||||
|
++Fout; |
||||||
|
} |
||||||
|
} |
||||||
|
} else |
||||||
|
{ |
||||||
|
kiss_fft_cpx * Fout_beg = Fout; |
||||||
|
for (i=0;i<N;i++) |
||||||
|
{ |
||||||
|
Fout = Fout_beg + i*mm; |
||||||
|
tw3 = tw2 = tw1 = st->twiddles; |
||||||
|
for (j=0;j<m;j++) |
||||||
|
{ |
||||||
|
C_MUL4(scratch[0],Fout[m] , *tw1 ); |
||||||
|
C_MUL4(scratch[1],Fout[m2] , *tw2 ); |
||||||
|
C_MUL4(scratch[2],Fout[m3] , *tw3 ); |
||||||
|
|
||||||
|
Fout->r = PSHR16(Fout->r, 2); |
||||||
|
Fout->i = PSHR16(Fout->i, 2); |
||||||
|
C_SUB( scratch[5] , *Fout, scratch[1] ); |
||||||
|
C_ADDTO(*Fout, scratch[1]); |
||||||
|
C_ADD( scratch[3] , scratch[0] , scratch[2] ); |
||||||
|
C_SUB( scratch[4] , scratch[0] , scratch[2] ); |
||||||
|
Fout[m2].r = PSHR16(Fout[m2].r, 2); |
||||||
|
Fout[m2].i = PSHR16(Fout[m2].i, 2); |
||||||
|
C_SUB( Fout[m2], *Fout, scratch[3] ); |
||||||
|
tw1 += fstride; |
||||||
|
tw2 += fstride*2; |
||||||
|
tw3 += fstride*3; |
||||||
|
C_ADDTO( *Fout , scratch[3] ); |
||||||
|
|
||||||
|
Fout[m].r = scratch[5].r + scratch[4].i; |
||||||
|
Fout[m].i = scratch[5].i - scratch[4].r; |
||||||
|
Fout[m3].r = scratch[5].r - scratch[4].i; |
||||||
|
Fout[m3].i = scratch[5].i + scratch[4].r; |
||||||
|
++Fout; |
||||||
|
} |
||||||
|
} |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
static void kf_bfly3( |
||||||
|
kiss_fft_cpx * Fout, |
||||||
|
const size_t fstride, |
||||||
|
const kiss_fft_cfg st, |
||||||
|
size_t m |
||||||
|
) |
||||||
|
{ |
||||||
|
size_t k=m; |
||||||
|
const size_t m2 = 2*m; |
||||||
|
kiss_fft_cpx *tw1,*tw2; |
||||||
|
kiss_fft_cpx scratch[5]; |
||||||
|
kiss_fft_cpx epi3; |
||||||
|
epi3 = st->twiddles[fstride*m]; |
||||||
|
|
||||||
|
tw1=tw2=st->twiddles; |
||||||
|
|
||||||
|
do{ |
||||||
|
if (!st->inverse) { |
||||||
|
C_FIXDIV(*Fout,3); C_FIXDIV(Fout[m],3); C_FIXDIV(Fout[m2],3); |
||||||
|
} |
||||||
|
|
||||||
|
C_MUL(scratch[1],Fout[m] , *tw1); |
||||||
|
C_MUL(scratch[2],Fout[m2] , *tw2); |
||||||
|
|
||||||
|
C_ADD(scratch[3],scratch[1],scratch[2]); |
||||||
|
C_SUB(scratch[0],scratch[1],scratch[2]); |
||||||
|
tw1 += fstride; |
||||||
|
tw2 += fstride*2; |
||||||
|
|
||||||
|
Fout[m].r = Fout->r - HALF_OF(scratch[3].r); |
||||||
|
Fout[m].i = Fout->i - HALF_OF(scratch[3].i); |
||||||
|
|
||||||
|
C_MULBYSCALAR( scratch[0] , epi3.i ); |
||||||
|
|
||||||
|
C_ADDTO(*Fout,scratch[3]); |
||||||
|
|
||||||
|
Fout[m2].r = Fout[m].r + scratch[0].i; |
||||||
|
Fout[m2].i = Fout[m].i - scratch[0].r; |
||||||
|
|
||||||
|
Fout[m].r -= scratch[0].i; |
||||||
|
Fout[m].i += scratch[0].r; |
||||||
|
|
||||||
|
++Fout; |
||||||
|
}while(--k); |
||||||
|
} |
||||||
|
|
||||||
|
static void kf_bfly5( |
||||||
|
kiss_fft_cpx * Fout, |
||||||
|
const size_t fstride, |
||||||
|
const kiss_fft_cfg st, |
||||||
|
int m |
||||||
|
) |
||||||
|
{ |
||||||
|
kiss_fft_cpx *Fout0,*Fout1,*Fout2,*Fout3,*Fout4; |
||||||
|
int u; |
||||||
|
kiss_fft_cpx scratch[13]; |
||||||
|
kiss_fft_cpx * twiddles = st->twiddles; |
||||||
|
kiss_fft_cpx *tw; |
||||||
|
kiss_fft_cpx ya,yb; |
||||||
|
ya = twiddles[fstride*m]; |
||||||
|
yb = twiddles[fstride*2*m]; |
||||||
|
|
||||||
|
Fout0=Fout; |
||||||
|
Fout1=Fout0+m; |
||||||
|
Fout2=Fout0+2*m; |
||||||
|
Fout3=Fout0+3*m; |
||||||
|
Fout4=Fout0+4*m; |
||||||
|
|
||||||
|
tw=st->twiddles; |
||||||
|
for ( u=0; u<m; ++u ) { |
||||||
|
if (!st->inverse) { |
||||||
|
C_FIXDIV( *Fout0,5); C_FIXDIV( *Fout1,5); C_FIXDIV( *Fout2,5); C_FIXDIV( *Fout3,5); C_FIXDIV( *Fout4,5); |
||||||
|
} |
||||||
|
scratch[0] = *Fout0; |
||||||
|
|
||||||
|
C_MUL(scratch[1] ,*Fout1, tw[u*fstride]); |
||||||
|
C_MUL(scratch[2] ,*Fout2, tw[2*u*fstride]); |
||||||
|
C_MUL(scratch[3] ,*Fout3, tw[3*u*fstride]); |
||||||
|
C_MUL(scratch[4] ,*Fout4, tw[4*u*fstride]); |
||||||
|
|
||||||
|
C_ADD( scratch[7],scratch[1],scratch[4]); |
||||||
|
C_SUB( scratch[10],scratch[1],scratch[4]); |
||||||
|
C_ADD( scratch[8],scratch[2],scratch[3]); |
||||||
|
C_SUB( scratch[9],scratch[2],scratch[3]); |
||||||
|
|
||||||
|
Fout0->r += scratch[7].r + scratch[8].r; |
||||||
|
Fout0->i += scratch[7].i + scratch[8].i; |
||||||
|
|
||||||
|
scratch[5].r = scratch[0].r + S_MUL(scratch[7].r,ya.r) + S_MUL(scratch[8].r,yb.r); |
||||||
|
scratch[5].i = scratch[0].i + S_MUL(scratch[7].i,ya.r) + S_MUL(scratch[8].i,yb.r); |
||||||
|
|
||||||
|
scratch[6].r = S_MUL(scratch[10].i,ya.i) + S_MUL(scratch[9].i,yb.i); |
||||||
|
scratch[6].i = -S_MUL(scratch[10].r,ya.i) - S_MUL(scratch[9].r,yb.i); |
||||||
|
|
||||||
|
C_SUB(*Fout1,scratch[5],scratch[6]); |
||||||
|
C_ADD(*Fout4,scratch[5],scratch[6]); |
||||||
|
|
||||||
|
scratch[11].r = scratch[0].r + S_MUL(scratch[7].r,yb.r) + S_MUL(scratch[8].r,ya.r); |
||||||
|
scratch[11].i = scratch[0].i + S_MUL(scratch[7].i,yb.r) + S_MUL(scratch[8].i,ya.r); |
||||||
|
scratch[12].r = - S_MUL(scratch[10].i,yb.i) + S_MUL(scratch[9].i,ya.i); |
||||||
|
scratch[12].i = S_MUL(scratch[10].r,yb.i) - S_MUL(scratch[9].r,ya.i); |
||||||
|
|
||||||
|
C_ADD(*Fout2,scratch[11],scratch[12]); |
||||||
|
C_SUB(*Fout3,scratch[11],scratch[12]); |
||||||
|
|
||||||
|
++Fout0;++Fout1;++Fout2;++Fout3;++Fout4; |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
/* perform the butterfly for one stage of a mixed radix FFT */ |
||||||
|
static void kf_bfly_generic( |
||||||
|
kiss_fft_cpx * Fout, |
||||||
|
const size_t fstride, |
||||||
|
const kiss_fft_cfg st, |
||||||
|
int m, |
||||||
|
int p |
||||||
|
) |
||||||
|
{ |
||||||
|
int u,k,q1,q; |
||||||
|
kiss_fft_cpx * twiddles = st->twiddles; |
||||||
|
kiss_fft_cpx t; |
||||||
|
kiss_fft_cpx scratchbuf[17]; |
||||||
|
int Norig = st->nfft; |
||||||
|
|
||||||
|
/*CHECKBUF(scratchbuf,nscratchbuf,p);*/ |
||||||
|
if (p>17) |
||||||
|
speex_fatal("KissFFT: max radix supported is 17"); |
||||||
|
|
||||||
|
for ( u=0; u<m; ++u ) { |
||||||
|
k=u; |
||||||
|
for ( q1=0 ; q1<p ; ++q1 ) { |
||||||
|
scratchbuf[q1] = Fout[ k ]; |
||||||
|
if (!st->inverse) { |
||||||
|
C_FIXDIV(scratchbuf[q1],p); |
||||||
|
} |
||||||
|
k += m; |
||||||
|
} |
||||||
|
|
||||||
|
k=u; |
||||||
|
for ( q1=0 ; q1<p ; ++q1 ) { |
||||||
|
int twidx=0; |
||||||
|
Fout[ k ] = scratchbuf[0]; |
||||||
|
for (q=1;q<p;++q ) { |
||||||
|
twidx += fstride * k; |
||||||
|
if (twidx>=Norig) twidx-=Norig; |
||||||
|
C_MUL(t,scratchbuf[q] , twiddles[twidx] ); |
||||||
|
C_ADDTO( Fout[ k ] ,t); |
||||||
|
} |
||||||
|
k += m; |
||||||
|
} |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
static |
||||||
|
void kf_shuffle( |
||||||
|
kiss_fft_cpx * Fout, |
||||||
|
const kiss_fft_cpx * f, |
||||||
|
const size_t fstride, |
||||||
|
int in_stride, |
||||||
|
int * factors, |
||||||
|
const kiss_fft_cfg st |
||||||
|
) |
||||||
|
{ |
||||||
|
const int p=*factors++; /* the radix */ |
||||||
|
const int m=*factors++; /* stage's fft length/p */ |
||||||
|
|
||||||
|
/*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/ |
||||||
|
if (m==1) |
||||||
|
{ |
||||||
|
int j; |
||||||
|
for (j=0;j<p;j++) |
||||||
|
{ |
||||||
|
Fout[j] = *f; |
||||||
|
f += fstride*in_stride; |
||||||
|
} |
||||||
|
} else { |
||||||
|
int j; |
||||||
|
for (j=0;j<p;j++) |
||||||
|
{ |
||||||
|
kf_shuffle( Fout , f, fstride*p, in_stride, factors,st); |
||||||
|
f += fstride*in_stride; |
||||||
|
Fout += m; |
||||||
|
} |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
static |
||||||
|
void kf_work( |
||||||
|
kiss_fft_cpx * Fout, |
||||||
|
const kiss_fft_cpx * f, |
||||||
|
const size_t fstride, |
||||||
|
int in_stride, |
||||||
|
int * factors, |
||||||
|
const kiss_fft_cfg st, |
||||||
|
int N, |
||||||
|
int s2, |
||||||
|
int m2 |
||||||
|
) |
||||||
|
{ |
||||||
|
int i; |
||||||
|
kiss_fft_cpx * Fout_beg=Fout; |
||||||
|
const int p=*factors++; /* the radix */ |
||||||
|
const int m=*factors++; /* stage's fft length/p */ |
||||||
|
#if 0 |
||||||
|
/*printf ("fft %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N);*/ |
||||||
|
if (m==1) |
||||||
|
{ |
||||||
|
/* int j;
|
||||||
|
for (j=0;j<p;j++) |
||||||
|
{ |
||||||
|
Fout[j] = *f; |
||||||
|
f += fstride*in_stride; |
||||||
|
}*/ |
||||||
|
} else { |
||||||
|
int j; |
||||||
|
for (j=0;j<p;j++) |
||||||
|
{ |
||||||
|
kf_work( Fout , f, fstride*p, in_stride, factors,st, N*p, fstride*in_stride, m); |
||||||
|
f += fstride*in_stride; |
||||||
|
Fout += m; |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
Fout=Fout_beg; |
||||||
|
|
||||||
|
switch (p) { |
||||||
|
case 2: kf_bfly2(Fout,fstride,st,m); break; |
||||||
|
case 3: kf_bfly3(Fout,fstride,st,m); break; |
||||||
|
case 4: kf_bfly4(Fout,fstride,st,m); break; |
||||||
|
case 5: kf_bfly5(Fout,fstride,st,m); break; |
||||||
|
default: kf_bfly_generic(Fout,fstride,st,m,p); break; |
||||||
|
} |
||||||
|
#else |
||||||
|
/*printf ("fft %d %d %d %d %d %d %d\n", p*m, m, p, s2, fstride*in_stride, N, m2);*/ |
||||||
|
if (m==1) |
||||||
|
{ |
||||||
|
/*for (i=0;i<N;i++)
|
||||||
|
{ |
||||||
|
int j; |
||||||
|
Fout = Fout_beg+i*m2; |
||||||
|
const kiss_fft_cpx * f2 = f+i*s2; |
||||||
|
for (j=0;j<p;j++) |
||||||
|
{ |
||||||
|
*Fout++ = *f2; |
||||||
|
f2 += fstride*in_stride; |
||||||
|
} |
||||||
|
}*/ |
||||||
|
}else{ |
||||||
|
kf_work( Fout , f, fstride*p, in_stride, factors,st, N*p, fstride*in_stride, m); |
||||||
|
} |
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
switch (p) { |
||||||
|
case 2: kf_bfly2(Fout,fstride,st,m, N, m2); break; |
||||||
|
case 3: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly3(Fout,fstride,st,m);} break; |
||||||
|
case 4: kf_bfly4(Fout,fstride,st,m, N, m2); break; |
||||||
|
case 5: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly5(Fout,fstride,st,m);} break; |
||||||
|
default: for (i=0;i<N;i++){Fout=Fout_beg+i*m2; kf_bfly_generic(Fout,fstride,st,m,p);} break; |
||||||
|
} |
||||||
|
#endif |
||||||
|
} |
||||||
|
|
||||||
|
/* facbuf is populated by p1,m1,p2,m2, ...
|
||||||
|
where |
||||||
|
p[i] * m[i] = m[i-1] |
||||||
|
m0 = n */ |
||||||
|
static |
||||||
|
void kf_factor(int n,int * facbuf) |
||||||
|
{ |
||||||
|
int p=4; |
||||||
|
|
||||||
|
/*factor out powers of 4, powers of 2, then any remaining primes */ |
||||||
|
do { |
||||||
|
while (n % p) { |
||||||
|
switch (p) { |
||||||
|
case 4: p = 2; break; |
||||||
|
case 2: p = 3; break; |
||||||
|
default: p += 2; break; |
||||||
|
} |
||||||
|
if (p>32000 || (spx_int32_t)p*(spx_int32_t)p > n) |
||||||
|
p = n; /* no more factors, skip to end */ |
||||||
|
} |
||||||
|
n /= p; |
||||||
|
*facbuf++ = p; |
||||||
|
*facbuf++ = n; |
||||||
|
} while (n > 1); |
||||||
|
} |
||||||
|
/*
|
||||||
|
* |
||||||
|
* User-callable function to allocate all necessary storage space for the fft. |
||||||
|
* |
||||||
|
* The return value is a contiguous block of memory, allocated with malloc. As such, |
||||||
|
* It can be freed with free(), rather than a kiss_fft-specific function. |
||||||
|
* */ |
||||||
|
kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem ) |
||||||
|
{ |
||||||
|
kiss_fft_cfg st=NULL; |
||||||
|
size_t memneeded = sizeof(struct kiss_fft_state) |
||||||
|
+ sizeof(kiss_fft_cpx)*(nfft-1); /* twiddle factors*/ |
||||||
|
|
||||||
|
if ( lenmem==NULL ) { |
||||||
|
st = ( kiss_fft_cfg)KISS_FFT_MALLOC( memneeded ); |
||||||
|
}else{ |
||||||
|
if (mem != NULL && *lenmem >= memneeded) |
||||||
|
st = (kiss_fft_cfg)mem; |
||||||
|
*lenmem = memneeded; |
||||||
|
} |
||||||
|
if (st) { |
||||||
|
int i; |
||||||
|
st->nfft=nfft; |
||||||
|
st->inverse = inverse_fft; |
||||||
|
#ifdef FIXED_POINT |
||||||
|
for (i=0;i<nfft;++i) { |
||||||
|
spx_word32_t phase = i; |
||||||
|
if (!st->inverse) |
||||||
|
phase = -phase; |
||||||
|
kf_cexp2(st->twiddles+i, DIV32(SHL32(phase,17),nfft)); |
||||||
|
} |
||||||
|
#else |
||||||
|
for (i=0;i<nfft;++i) { |
||||||
|
const double pi=3.14159265358979323846264338327; |
||||||
|
double phase = ( -2*pi /nfft ) * i; |
||||||
|
if (st->inverse) |
||||||
|
phase *= -1; |
||||||
|
kf_cexp(st->twiddles+i, phase ); |
||||||
|
} |
||||||
|
#endif |
||||||
|
kf_factor(nfft,st->factors); |
||||||
|
} |
||||||
|
return st; |
||||||
|
} |
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
void kiss_fft_stride(kiss_fft_cfg st,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int in_stride) |
||||||
|
{ |
||||||
|
if (fin == fout) |
||||||
|
{ |
||||||
|
speex_fatal("In-place FFT not supported"); |
||||||
|
/*CHECKBUF(tmpbuf,ntmpbuf,st->nfft);
|
||||||
|
kf_work(tmpbuf,fin,1,in_stride, st->factors,st); |
||||||
|
SPEEX_MOVE(fout,tmpbuf,st->nfft);*/ |
||||||
|
} else { |
||||||
|
kf_shuffle( fout, fin, 1,in_stride, st->factors,st); |
||||||
|
kf_work( fout, fin, 1,in_stride, st->factors,st, 1, in_stride, 1); |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout) |
||||||
|
{ |
||||||
|
kiss_fft_stride(cfg,fin,fout,1); |
||||||
|
} |
||||||
|
|
||||||
@ -0,0 +1,108 @@ |
|||||||
|
#ifndef KISS_FFT_H |
||||||
|
#define KISS_FFT_H |
||||||
|
|
||||||
|
#include <stdlib.h> |
||||||
|
#include <math.h> |
||||||
|
#include "arch.h" |
||||||
|
|
||||||
|
#ifdef __cplusplus |
||||||
|
extern "C" { |
||||||
|
#endif |
||||||
|
|
||||||
|
/*
|
||||||
|
ATTENTION! |
||||||
|
If you would like a : |
||||||
|
-- a utility that will handle the caching of fft objects |
||||||
|
-- real-only (no imaginary time component ) FFT |
||||||
|
-- a multi-dimensional FFT |
||||||
|
-- a command-line utility to perform ffts |
||||||
|
-- a command-line utility to perform fast-convolution filtering |
||||||
|
|
||||||
|
Then see kfc.h kiss_fftr.h kiss_fftnd.h fftutil.c kiss_fastfir.c |
||||||
|
in the tools/ directory. |
||||||
|
*/ |
||||||
|
|
||||||
|
#ifdef USE_SIMD |
||||||
|
# include <xmmintrin.h> |
||||||
|
# define kiss_fft_scalar __m128 |
||||||
|
#define KISS_FFT_MALLOC(nbytes) memalign(16,nbytes) |
||||||
|
#else |
||||||
|
#define KISS_FFT_MALLOC speex_alloc |
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
#ifdef FIXED_POINT |
||||||
|
#include "arch.h" |
||||||
|
# define kiss_fft_scalar spx_int16_t |
||||||
|
#else |
||||||
|
# ifndef kiss_fft_scalar |
||||||
|
/* default is float */ |
||||||
|
# define kiss_fft_scalar float |
||||||
|
# endif |
||||||
|
#endif |
||||||
|
|
||||||
|
typedef struct { |
||||||
|
kiss_fft_scalar r; |
||||||
|
kiss_fft_scalar i; |
||||||
|
}kiss_fft_cpx; |
||||||
|
|
||||||
|
typedef struct kiss_fft_state* kiss_fft_cfg; |
||||||
|
|
||||||
|
/*
|
||||||
|
* kiss_fft_alloc |
||||||
|
* |
||||||
|
* Initialize a FFT (or IFFT) algorithm's cfg/state buffer. |
||||||
|
* |
||||||
|
* typical usage: kiss_fft_cfg mycfg=kiss_fft_alloc(1024,0,NULL,NULL); |
||||||
|
* |
||||||
|
* The return value from fft_alloc is a cfg buffer used internally |
||||||
|
* by the fft routine or NULL. |
||||||
|
* |
||||||
|
* If lenmem is NULL, then kiss_fft_alloc will allocate a cfg buffer using malloc. |
||||||
|
* The returned value should be free()d when done to avoid memory leaks. |
||||||
|
* |
||||||
|
* The state can be placed in a user supplied buffer 'mem': |
||||||
|
* If lenmem is not NULL and mem is not NULL and *lenmem is large enough, |
||||||
|
* then the function places the cfg in mem and the size used in *lenmem |
||||||
|
* and returns mem. |
||||||
|
* |
||||||
|
* If lenmem is not NULL and ( mem is NULL or *lenmem is not large enough), |
||||||
|
* then the function returns NULL and places the minimum cfg |
||||||
|
* buffer size in *lenmem. |
||||||
|
* */ |
||||||
|
|
||||||
|
kiss_fft_cfg kiss_fft_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem); |
||||||
|
|
||||||
|
/*
|
||||||
|
* kiss_fft(cfg,in_out_buf) |
||||||
|
* |
||||||
|
* Perform an FFT on a complex input buffer. |
||||||
|
* for a forward FFT, |
||||||
|
* fin should be f[0] , f[1] , ... ,f[nfft-1] |
||||||
|
* fout will be F[0] , F[1] , ... ,F[nfft-1] |
||||||
|
* Note that each element is complex and can be accessed like |
||||||
|
f[k].r and f[k].i |
||||||
|
* */ |
||||||
|
void kiss_fft(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout); |
||||||
|
|
||||||
|
/*
|
||||||
|
A more generic version of the above function. It reads its input from every Nth sample. |
||||||
|
* */ |
||||||
|
void kiss_fft_stride(kiss_fft_cfg cfg,const kiss_fft_cpx *fin,kiss_fft_cpx *fout,int fin_stride); |
||||||
|
|
||||||
|
/* If kiss_fft_alloc allocated a buffer, it is one contiguous
|
||||||
|
buffer and can be simply free()d when no longer needed*/ |
||||||
|
#define kiss_fft_free speex_free |
||||||
|
|
||||||
|
/*
|
||||||
|
Cleans up some memory that gets managed internally. Not necessary to call, but it might clean up |
||||||
|
your compiler output to call this before you exit. |
||||||
|
*/ |
||||||
|
void kiss_fft_cleanup(void); |
||||||
|
|
||||||
|
|
||||||
|
#ifdef __cplusplus |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
#endif |
||||||
@ -0,0 +1,297 @@ |
|||||||
|
/*
|
||||||
|
Copyright (c) 2003-2004, Mark Borgerding |
||||||
|
|
||||||
|
All rights reserved. |
||||||
|
|
||||||
|
Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: |
||||||
|
|
||||||
|
* Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. |
||||||
|
* Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. |
||||||
|
* Neither the author nor the names of any contributors may be used to endorse or promote products derived from this software without specific prior written permission. |
||||||
|
|
||||||
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
||||||
|
*/ |
||||||
|
|
||||||
|
#ifdef HAVE_CONFIG_H |
||||||
|
#include "config.h" |
||||||
|
#endif |
||||||
|
|
||||||
|
#include "os_support.h" |
||||||
|
#include "kiss_fftr.h" |
||||||
|
#include "_kiss_fft_guts.h" |
||||||
|
|
||||||
|
struct kiss_fftr_state{ |
||||||
|
kiss_fft_cfg substate; |
||||||
|
kiss_fft_cpx * tmpbuf; |
||||||
|
kiss_fft_cpx * super_twiddles; |
||||||
|
#ifdef USE_SIMD |
||||||
|
long pad; |
||||||
|
#endif |
||||||
|
}; |
||||||
|
|
||||||
|
kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem,size_t * lenmem) |
||||||
|
{ |
||||||
|
int i; |
||||||
|
kiss_fftr_cfg st = NULL; |
||||||
|
size_t subsize, memneeded; |
||||||
|
|
||||||
|
if (nfft & 1) { |
||||||
|
speex_warning("Real FFT optimization must be even.\n"); |
||||||
|
return NULL; |
||||||
|
} |
||||||
|
nfft >>= 1; |
||||||
|
|
||||||
|
kiss_fft_alloc (nfft, inverse_fft, NULL, &subsize); |
||||||
|
memneeded = sizeof(struct kiss_fftr_state) + subsize + sizeof(kiss_fft_cpx) * ( nfft * 2); |
||||||
|
|
||||||
|
if (lenmem == NULL) { |
||||||
|
st = (kiss_fftr_cfg) KISS_FFT_MALLOC (memneeded); |
||||||
|
} else { |
||||||
|
if (*lenmem >= memneeded) |
||||||
|
st = (kiss_fftr_cfg) mem; |
||||||
|
*lenmem = memneeded; |
||||||
|
} |
||||||
|
if (!st) |
||||||
|
return NULL; |
||||||
|
|
||||||
|
st->substate = (kiss_fft_cfg) (st + 1); /*just beyond kiss_fftr_state struct */ |
||||||
|
st->tmpbuf = (kiss_fft_cpx *) (((char *) st->substate) + subsize); |
||||||
|
st->super_twiddles = st->tmpbuf + nfft; |
||||||
|
kiss_fft_alloc(nfft, inverse_fft, st->substate, &subsize); |
||||||
|
|
||||||
|
#ifdef FIXED_POINT |
||||||
|
for (i=0;i<nfft;++i) { |
||||||
|
spx_word32_t phase = i+(nfft>>1); |
||||||
|
if (!inverse_fft) |
||||||
|
phase = -phase; |
||||||
|
kf_cexp2(st->super_twiddles+i, DIV32(SHL32(phase,16),nfft)); |
||||||
|
} |
||||||
|
#else |
||||||
|
for (i=0;i<nfft;++i) { |
||||||
|
const double pi=3.14159265358979323846264338327; |
||||||
|
double phase = pi*(((double)i) /nfft + .5); |
||||||
|
if (!inverse_fft) |
||||||
|
phase = -phase; |
||||||
|
kf_cexp(st->super_twiddles+i, phase ); |
||||||
|
} |
||||||
|
#endif |
||||||
|
return st; |
||||||
|
} |
||||||
|
|
||||||
|
void kiss_fftr(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata) |
||||||
|
{ |
||||||
|
/* input buffer timedata is stored row-wise */ |
||||||
|
int k,ncfft; |
||||||
|
kiss_fft_cpx fpnk,fpk,f1k,f2k,tw,tdc; |
||||||
|
|
||||||
|
if ( st->substate->inverse) { |
||||||
|
speex_fatal("kiss fft usage error: improper alloc\n"); |
||||||
|
} |
||||||
|
|
||||||
|
ncfft = st->substate->nfft; |
||||||
|
|
||||||
|
/*perform the parallel fft of two real signals packed in real,imag*/ |
||||||
|
kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf ); |
||||||
|
/* The real part of the DC element of the frequency spectrum in st->tmpbuf
|
||||||
|
* contains the sum of the even-numbered elements of the input time sequence |
||||||
|
* The imag part is the sum of the odd-numbered elements |
||||||
|
* |
||||||
|
* The sum of tdc.r and tdc.i is the sum of the input time sequence. |
||||||
|
* yielding DC of input time sequence |
||||||
|
* The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1... |
||||||
|
* yielding Nyquist bin of input time sequence |
||||||
|
*/ |
||||||
|
|
||||||
|
tdc.r = st->tmpbuf[0].r; |
||||||
|
tdc.i = st->tmpbuf[0].i; |
||||||
|
C_FIXDIV(tdc,2); |
||||||
|
CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i); |
||||||
|
CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i); |
||||||
|
freqdata[0].r = tdc.r + tdc.i; |
||||||
|
freqdata[ncfft].r = tdc.r - tdc.i; |
||||||
|
#ifdef USE_SIMD |
||||||
|
freqdata[ncfft].i = freqdata[0].i = _mm_set1_ps(0); |
||||||
|
#else |
||||||
|
freqdata[ncfft].i = freqdata[0].i = 0; |
||||||
|
#endif |
||||||
|
|
||||||
|
for ( k=1;k <= ncfft/2 ; ++k ) { |
||||||
|
fpk = st->tmpbuf[k]; |
||||||
|
fpnk.r = st->tmpbuf[ncfft-k].r; |
||||||
|
fpnk.i = - st->tmpbuf[ncfft-k].i; |
||||||
|
C_FIXDIV(fpk,2); |
||||||
|
C_FIXDIV(fpnk,2); |
||||||
|
|
||||||
|
C_ADD( f1k, fpk , fpnk ); |
||||||
|
C_SUB( f2k, fpk , fpnk ); |
||||||
|
C_MUL( tw , f2k , st->super_twiddles[k]); |
||||||
|
|
||||||
|
freqdata[k].r = HALF_OF(f1k.r + tw.r); |
||||||
|
freqdata[k].i = HALF_OF(f1k.i + tw.i); |
||||||
|
freqdata[ncfft-k].r = HALF_OF(f1k.r - tw.r); |
||||||
|
freqdata[ncfft-k].i = HALF_OF(tw.i - f1k.i); |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
void kiss_fftri(kiss_fftr_cfg st,const kiss_fft_cpx *freqdata, kiss_fft_scalar *timedata) |
||||||
|
{ |
||||||
|
/* input buffer timedata is stored row-wise */ |
||||||
|
int k, ncfft; |
||||||
|
|
||||||
|
if (st->substate->inverse == 0) { |
||||||
|
speex_fatal("kiss fft usage error: improper alloc\n"); |
||||||
|
} |
||||||
|
|
||||||
|
ncfft = st->substate->nfft; |
||||||
|
|
||||||
|
st->tmpbuf[0].r = freqdata[0].r + freqdata[ncfft].r; |
||||||
|
st->tmpbuf[0].i = freqdata[0].r - freqdata[ncfft].r; |
||||||
|
/*C_FIXDIV(st->tmpbuf[0],2);*/ |
||||||
|
|
||||||
|
for (k = 1; k <= ncfft / 2; ++k) { |
||||||
|
kiss_fft_cpx fk, fnkc, fek, fok, tmp; |
||||||
|
fk = freqdata[k]; |
||||||
|
fnkc.r = freqdata[ncfft - k].r; |
||||||
|
fnkc.i = -freqdata[ncfft - k].i; |
||||||
|
/*C_FIXDIV( fk , 2 );
|
||||||
|
C_FIXDIV( fnkc , 2 );*/ |
||||||
|
|
||||||
|
C_ADD (fek, fk, fnkc); |
||||||
|
C_SUB (tmp, fk, fnkc); |
||||||
|
C_MUL (fok, tmp, st->super_twiddles[k]); |
||||||
|
C_ADD (st->tmpbuf[k], fek, fok); |
||||||
|
C_SUB (st->tmpbuf[ncfft - k], fek, fok); |
||||||
|
#ifdef USE_SIMD |
||||||
|
st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0); |
||||||
|
#else |
||||||
|
st->tmpbuf[ncfft - k].i *= -1; |
||||||
|
#endif |
||||||
|
} |
||||||
|
kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata); |
||||||
|
} |
||||||
|
|
||||||
|
void kiss_fftr2(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_scalar *freqdata) |
||||||
|
{ |
||||||
|
/* input buffer timedata is stored row-wise */ |
||||||
|
int k,ncfft; |
||||||
|
kiss_fft_cpx f2k,tdc; |
||||||
|
spx_word32_t f1kr, f1ki, twr, twi; |
||||||
|
|
||||||
|
if ( st->substate->inverse) { |
||||||
|
speex_fatal("kiss fft usage error: improper alloc\n"); |
||||||
|
} |
||||||
|
|
||||||
|
ncfft = st->substate->nfft; |
||||||
|
|
||||||
|
/*perform the parallel fft of two real signals packed in real,imag*/ |
||||||
|
kiss_fft( st->substate , (const kiss_fft_cpx*)timedata, st->tmpbuf ); |
||||||
|
/* The real part of the DC element of the frequency spectrum in st->tmpbuf
|
||||||
|
* contains the sum of the even-numbered elements of the input time sequence |
||||||
|
* The imag part is the sum of the odd-numbered elements |
||||||
|
* |
||||||
|
* The sum of tdc.r and tdc.i is the sum of the input time sequence. |
||||||
|
* yielding DC of input time sequence |
||||||
|
* The difference of tdc.r - tdc.i is the sum of the input (dot product) [1,-1,1,-1... |
||||||
|
* yielding Nyquist bin of input time sequence |
||||||
|
*/ |
||||||
|
|
||||||
|
tdc.r = st->tmpbuf[0].r; |
||||||
|
tdc.i = st->tmpbuf[0].i; |
||||||
|
C_FIXDIV(tdc,2); |
||||||
|
CHECK_OVERFLOW_OP(tdc.r ,+, tdc.i); |
||||||
|
CHECK_OVERFLOW_OP(tdc.r ,-, tdc.i); |
||||||
|
freqdata[0] = tdc.r + tdc.i; |
||||||
|
freqdata[2*ncfft-1] = tdc.r - tdc.i; |
||||||
|
|
||||||
|
for ( k=1;k <= ncfft/2 ; ++k ) |
||||||
|
{ |
||||||
|
/*fpk = st->tmpbuf[k];
|
||||||
|
fpnk.r = st->tmpbuf[ncfft-k].r; |
||||||
|
fpnk.i = - st->tmpbuf[ncfft-k].i; |
||||||
|
C_FIXDIV(fpk,2); |
||||||
|
C_FIXDIV(fpnk,2); |
||||||
|
|
||||||
|
C_ADD( f1k, fpk , fpnk ); |
||||||
|
C_SUB( f2k, fpk , fpnk ); |
||||||
|
|
||||||
|
C_MUL( tw , f2k , st->super_twiddles[k]); |
||||||
|
|
||||||
|
freqdata[2*k-1] = HALF_OF(f1k.r + tw.r); |
||||||
|
freqdata[2*k] = HALF_OF(f1k.i + tw.i); |
||||||
|
freqdata[2*(ncfft-k)-1] = HALF_OF(f1k.r - tw.r); |
||||||
|
freqdata[2*(ncfft-k)] = HALF_OF(tw.i - f1k.i); |
||||||
|
*/ |
||||||
|
|
||||||
|
/*f1k.r = PSHR32(ADD32(EXTEND32(st->tmpbuf[k].r), EXTEND32(st->tmpbuf[ncfft-k].r)),1);
|
||||||
|
f1k.i = PSHR32(SUB32(EXTEND32(st->tmpbuf[k].i), EXTEND32(st->tmpbuf[ncfft-k].i)),1); |
||||||
|
f2k.r = PSHR32(SUB32(EXTEND32(st->tmpbuf[k].r), EXTEND32(st->tmpbuf[ncfft-k].r)),1); |
||||||
|
f2k.i = SHR32(ADD32(EXTEND32(st->tmpbuf[k].i), EXTEND32(st->tmpbuf[ncfft-k].i)),1); |
||||||
|
|
||||||
|
C_MUL( tw , f2k , st->super_twiddles[k]); |
||||||
|
|
||||||
|
freqdata[2*k-1] = HALF_OF(f1k.r + tw.r); |
||||||
|
freqdata[2*k] = HALF_OF(f1k.i + tw.i); |
||||||
|
freqdata[2*(ncfft-k)-1] = HALF_OF(f1k.r - tw.r); |
||||||
|
freqdata[2*(ncfft-k)] = HALF_OF(tw.i - f1k.i); |
||||||
|
*/ |
||||||
|
f2k.r = SHR32(SUB32(EXTEND32(st->tmpbuf[k].r), EXTEND32(st->tmpbuf[ncfft-k].r)),1); |
||||||
|
f2k.i = PSHR32(ADD32(EXTEND32(st->tmpbuf[k].i), EXTEND32(st->tmpbuf[ncfft-k].i)),1); |
||||||
|
|
||||||
|
f1kr = SHL32(ADD32(EXTEND32(st->tmpbuf[k].r), EXTEND32(st->tmpbuf[ncfft-k].r)),13); |
||||||
|
f1ki = SHL32(SUB32(EXTEND32(st->tmpbuf[k].i), EXTEND32(st->tmpbuf[ncfft-k].i)),13); |
||||||
|
|
||||||
|
twr = SHR32(SUB32(MULT16_16(f2k.r,st->super_twiddles[k].r),MULT16_16(f2k.i,st->super_twiddles[k].i)), 1); |
||||||
|
twi = SHR32(ADD32(MULT16_16(f2k.i,st->super_twiddles[k].r),MULT16_16(f2k.r,st->super_twiddles[k].i)), 1); |
||||||
|
|
||||||
|
#ifdef FIXED_POINT |
||||||
|
freqdata[2*k-1] = PSHR32(f1kr + twr, 15); |
||||||
|
freqdata[2*k] = PSHR32(f1ki + twi, 15); |
||||||
|
freqdata[2*(ncfft-k)-1] = PSHR32(f1kr - twr, 15); |
||||||
|
freqdata[2*(ncfft-k)] = PSHR32(twi - f1ki, 15); |
||||||
|
#else |
||||||
|
freqdata[2*k-1] = .5f*(f1kr + twr); |
||||||
|
freqdata[2*k] = .5f*(f1ki + twi); |
||||||
|
freqdata[2*(ncfft-k)-1] = .5f*(f1kr - twr); |
||||||
|
freqdata[2*(ncfft-k)] = .5f*(twi - f1ki); |
||||||
|
|
||||||
|
#endif |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
void kiss_fftri2(kiss_fftr_cfg st,const kiss_fft_scalar *freqdata,kiss_fft_scalar *timedata) |
||||||
|
{ |
||||||
|
/* input buffer timedata is stored row-wise */ |
||||||
|
int k, ncfft; |
||||||
|
|
||||||
|
if (st->substate->inverse == 0) { |
||||||
|
speex_fatal ("kiss fft usage error: improper alloc\n"); |
||||||
|
} |
||||||
|
|
||||||
|
ncfft = st->substate->nfft; |
||||||
|
|
||||||
|
st->tmpbuf[0].r = freqdata[0] + freqdata[2*ncfft-1]; |
||||||
|
st->tmpbuf[0].i = freqdata[0] - freqdata[2*ncfft-1]; |
||||||
|
/*C_FIXDIV(st->tmpbuf[0],2);*/ |
||||||
|
|
||||||
|
for (k = 1; k <= ncfft / 2; ++k) { |
||||||
|
kiss_fft_cpx fk, fnkc, fek, fok, tmp; |
||||||
|
fk.r = freqdata[2*k-1]; |
||||||
|
fk.i = freqdata[2*k]; |
||||||
|
fnkc.r = freqdata[2*(ncfft - k)-1]; |
||||||
|
fnkc.i = -freqdata[2*(ncfft - k)]; |
||||||
|
/*C_FIXDIV( fk , 2 );
|
||||||
|
C_FIXDIV( fnkc , 2 );*/ |
||||||
|
|
||||||
|
C_ADD (fek, fk, fnkc); |
||||||
|
C_SUB (tmp, fk, fnkc); |
||||||
|
C_MUL (fok, tmp, st->super_twiddles[k]); |
||||||
|
C_ADD (st->tmpbuf[k], fek, fok); |
||||||
|
C_SUB (st->tmpbuf[ncfft - k], fek, fok); |
||||||
|
#ifdef USE_SIMD |
||||||
|
st->tmpbuf[ncfft - k].i *= _mm_set1_ps(-1.0); |
||||||
|
#else |
||||||
|
st->tmpbuf[ncfft - k].i *= -1; |
||||||
|
#endif |
||||||
|
} |
||||||
|
kiss_fft (st->substate, st->tmpbuf, (kiss_fft_cpx *) timedata); |
||||||
|
} |
||||||
@ -0,0 +1,51 @@ |
|||||||
|
#ifndef KISS_FTR_H |
||||||
|
#define KISS_FTR_H |
||||||
|
|
||||||
|
#include "kiss_fft.h" |
||||||
|
#ifdef __cplusplus |
||||||
|
extern "C" { |
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
/*
|
||||||
|
|
||||||
|
Real optimized version can save about 45% cpu time vs. complex fft of a real seq. |
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
*/ |
||||||
|
|
||||||
|
typedef struct kiss_fftr_state *kiss_fftr_cfg; |
||||||
|
|
||||||
|
|
||||||
|
kiss_fftr_cfg kiss_fftr_alloc(int nfft,int inverse_fft,void * mem, size_t * lenmem); |
||||||
|
/*
|
||||||
|
nfft must be even |
||||||
|
|
||||||
|
If you don't care to allocate space, use mem = lenmem = NULL |
||||||
|
*/ |
||||||
|
|
||||||
|
|
||||||
|
void kiss_fftr(kiss_fftr_cfg cfg,const kiss_fft_scalar *timedata,kiss_fft_cpx *freqdata); |
||||||
|
/*
|
||||||
|
input timedata has nfft scalar points |
||||||
|
output freqdata has nfft/2+1 complex points |
||||||
|
*/ |
||||||
|
|
||||||
|
void kiss_fftr2(kiss_fftr_cfg st,const kiss_fft_scalar *timedata,kiss_fft_scalar *freqdata); |
||||||
|
|
||||||
|
void kiss_fftri(kiss_fftr_cfg cfg,const kiss_fft_cpx *freqdata,kiss_fft_scalar *timedata); |
||||||
|
|
||||||
|
void kiss_fftri2(kiss_fftr_cfg st,const kiss_fft_scalar *freqdata, kiss_fft_scalar *timedata); |
||||||
|
|
||||||
|
/*
|
||||||
|
input freqdata has nfft/2+1 complex points |
||||||
|
output timedata has nfft scalar points |
||||||
|
*/ |
||||||
|
|
||||||
|
#define kiss_fftr_free speex_free |
||||||
|
|
||||||
|
#ifdef __cplusplus |
||||||
|
} |
||||||
|
#endif |
||||||
|
#endif |
||||||
@ -0,0 +1,332 @@ |
|||||||
|
/* Copyright (C) 2002 Jean-Marc Valin */ |
||||||
|
/**
|
||||||
|
@file math_approx.h |
||||||
|
@brief Various math approximation functions for Speex |
||||||
|
*/ |
||||||
|
/*
|
||||||
|
Redistribution and use in source and binary forms, with or without |
||||||
|
modification, are permitted provided that the following conditions |
||||||
|
are met: |
||||||
|
|
||||||
|
- Redistributions of source code must retain the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer. |
||||||
|
|
||||||
|
- Redistributions in binary form must reproduce the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer in the |
||||||
|
documentation and/or other materials provided with the distribution. |
||||||
|
|
||||||
|
- Neither the name of the Xiph.org Foundation nor the names of its |
||||||
|
contributors may be used to endorse or promote products derived from |
||||||
|
this software without specific prior written permission. |
||||||
|
|
||||||
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
||||||
|
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
||||||
|
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
||||||
|
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR |
||||||
|
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
||||||
|
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
||||||
|
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
||||||
|
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
||||||
|
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
||||||
|
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
||||||
|
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
||||||
|
*/ |
||||||
|
|
||||||
|
#ifndef MATH_APPROX_H |
||||||
|
#define MATH_APPROX_H |
||||||
|
|
||||||
|
#include "arch.h" |
||||||
|
|
||||||
|
#ifndef FIXED_POINT |
||||||
|
|
||||||
|
#define spx_sqrt sqrt |
||||||
|
#define spx_acos acos |
||||||
|
#define spx_exp exp |
||||||
|
#define spx_cos_norm(x) (cos((.5f*M_PI)*(x))) |
||||||
|
#define spx_atan atan |
||||||
|
|
||||||
|
/** Generate a pseudo-random number */ |
||||||
|
static inline spx_word16_t speex_rand(spx_word16_t std, spx_int32_t *seed) |
||||||
|
{ |
||||||
|
const unsigned int jflone = 0x3f800000; |
||||||
|
const unsigned int jflmsk = 0x007fffff; |
||||||
|
union {int i; float f;} ran; |
||||||
|
*seed = 1664525 * *seed + 1013904223; |
||||||
|
ran.i = jflone | (jflmsk & *seed); |
||||||
|
ran.f -= 1.5; |
||||||
|
return 3.4642*std*ran.f; |
||||||
|
} |
||||||
|
|
||||||
|
|
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
static inline spx_int16_t spx_ilog2(spx_uint32_t x) |
||||||
|
{ |
||||||
|
int r=0; |
||||||
|
if (x>=(spx_int32_t)65536) |
||||||
|
{ |
||||||
|
x >>= 16; |
||||||
|
r += 16; |
||||||
|
} |
||||||
|
if (x>=256) |
||||||
|
{ |
||||||
|
x >>= 8; |
||||||
|
r += 8; |
||||||
|
} |
||||||
|
if (x>=16) |
||||||
|
{ |
||||||
|
x >>= 4; |
||||||
|
r += 4; |
||||||
|
} |
||||||
|
if (x>=4) |
||||||
|
{ |
||||||
|
x >>= 2; |
||||||
|
r += 2; |
||||||
|
} |
||||||
|
if (x>=2) |
||||||
|
{ |
||||||
|
r += 1; |
||||||
|
} |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
static inline spx_int16_t spx_ilog4(spx_uint32_t x) |
||||||
|
{ |
||||||
|
int r=0; |
||||||
|
if (x>=(spx_int32_t)65536) |
||||||
|
{ |
||||||
|
x >>= 16; |
||||||
|
r += 8; |
||||||
|
} |
||||||
|
if (x>=256) |
||||||
|
{ |
||||||
|
x >>= 8; |
||||||
|
r += 4; |
||||||
|
} |
||||||
|
if (x>=16) |
||||||
|
{ |
||||||
|
x >>= 4; |
||||||
|
r += 2; |
||||||
|
} |
||||||
|
if (x>=4) |
||||||
|
{ |
||||||
|
r += 1; |
||||||
|
} |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
#ifdef FIXED_POINT |
||||||
|
|
||||||
|
/** Generate a pseudo-random number */ |
||||||
|
static inline spx_word16_t speex_rand(spx_word16_t std, spx_int32_t *seed) |
||||||
|
{ |
||||||
|
spx_word32_t res; |
||||||
|
*seed = 1664525 * *seed + 1013904223; |
||||||
|
res = MULT16_16(EXTRACT16(SHR32(*seed,16)),std); |
||||||
|
return EXTRACT16(PSHR32(SUB32(res, SHR32(res, 3)),14)); |
||||||
|
} |
||||||
|
|
||||||
|
/* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25723*x^3 (for .25 < x < 1) */ |
||||||
|
/*#define C0 3634
|
||||||
|
#define C1 21173 |
||||||
|
#define C2 -12627 |
||||||
|
#define C3 4215*/ |
||||||
|
|
||||||
|
/* sqrt(x) ~= 0.22178 + 1.29227*x - 0.77070*x^2 + 0.25659*x^3 (for .25 < x < 1) */ |
||||||
|
#define C0 3634 |
||||||
|
#define C1 21173 |
||||||
|
#define C2 -12627 |
||||||
|
#define C3 4204 |
||||||
|
|
||||||
|
static inline spx_word16_t spx_sqrt(spx_word32_t x) |
||||||
|
{ |
||||||
|
int k; |
||||||
|
spx_word32_t rt; |
||||||
|
k = spx_ilog4(x)-6; |
||||||
|
x = VSHR32(x, (k<<1)); |
||||||
|
rt = ADD16(C0, MULT16_16_Q14(x, ADD16(C1, MULT16_16_Q14(x, ADD16(C2, MULT16_16_Q14(x, (C3))))))); |
||||||
|
rt = VSHR32(rt,7-k); |
||||||
|
return rt; |
||||||
|
} |
||||||
|
|
||||||
|
/* log(x) ~= -2.18151 + 4.20592*x - 2.88938*x^2 + 0.86535*x^3 (for .5 < x < 1) */ |
||||||
|
|
||||||
|
|
||||||
|
#define A1 16469 |
||||||
|
#define A2 2242 |
||||||
|
#define A3 1486 |
||||||
|
|
||||||
|
static inline spx_word16_t spx_acos(spx_word16_t x) |
||||||
|
{ |
||||||
|
int s=0; |
||||||
|
spx_word16_t ret; |
||||||
|
spx_word16_t sq; |
||||||
|
if (x<0) |
||||||
|
{ |
||||||
|
s=1; |
||||||
|
x = NEG16(x); |
||||||
|
} |
||||||
|
x = SUB16(16384,x); |
||||||
|
|
||||||
|
x = x >> 1; |
||||||
|
sq = MULT16_16_Q13(x, ADD16(A1, MULT16_16_Q13(x, ADD16(A2, MULT16_16_Q13(x, (A3)))))); |
||||||
|
ret = spx_sqrt(SHL32(EXTEND32(sq),13)); |
||||||
|
|
||||||
|
/*ret = spx_sqrt(67108864*(-1.6129e-04 + 2.0104e+00*f + 2.7373e-01*f*f + 1.8136e-01*f*f*f));*/ |
||||||
|
if (s) |
||||||
|
ret = SUB16(25736,ret); |
||||||
|
return ret; |
||||||
|
} |
||||||
|
|
||||||
|
|
||||||
|
#define K1 8192 |
||||||
|
#define K2 -4096 |
||||||
|
#define K3 340 |
||||||
|
#define K4 -10 |
||||||
|
|
||||||
|
static inline spx_word16_t spx_cos(spx_word16_t x) |
||||||
|
{ |
||||||
|
spx_word16_t x2; |
||||||
|
|
||||||
|
if (x<12868) |
||||||
|
{ |
||||||
|
x2 = MULT16_16_P13(x,x); |
||||||
|
return ADD32(K1, MULT16_16_P13(x2, ADD32(K2, MULT16_16_P13(x2, ADD32(K3, MULT16_16_P13(K4, x2)))))); |
||||||
|
} else { |
||||||
|
x = SUB16(25736,x); |
||||||
|
x2 = MULT16_16_P13(x,x); |
||||||
|
return SUB32(-K1, MULT16_16_P13(x2, ADD32(K2, MULT16_16_P13(x2, ADD32(K3, MULT16_16_P13(K4, x2)))))); |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
#define L1 32767 |
||||||
|
#define L2 -7651 |
||||||
|
#define L3 8277 |
||||||
|
#define L4 -626 |
||||||
|
|
||||||
|
static inline spx_word16_t _spx_cos_pi_2(spx_word16_t x) |
||||||
|
{ |
||||||
|
spx_word16_t x2; |
||||||
|
|
||||||
|
x2 = MULT16_16_P15(x,x); |
||||||
|
return ADD16(1,MIN16(32766,ADD32(SUB16(L1,x2), MULT16_16_P15(x2, ADD32(L2, MULT16_16_P15(x2, ADD32(L3, MULT16_16_P15(L4, x2)))))))); |
||||||
|
} |
||||||
|
|
||||||
|
static inline spx_word16_t spx_cos_norm(spx_word32_t x) |
||||||
|
{ |
||||||
|
x = x&0x0001ffff; |
||||||
|
if (x>SHL32(EXTEND32(1), 16)) |
||||||
|
x = SUB32(SHL32(EXTEND32(1), 17),x); |
||||||
|
if (x&0x00007fff) |
||||||
|
{ |
||||||
|
if (x<SHL32(EXTEND32(1), 15)) |
||||||
|
{ |
||||||
|
return _spx_cos_pi_2(EXTRACT16(x)); |
||||||
|
} else { |
||||||
|
return NEG32(_spx_cos_pi_2(EXTRACT16(65536-x))); |
||||||
|
} |
||||||
|
} else { |
||||||
|
if (x&0x0000ffff) |
||||||
|
return 0; |
||||||
|
else if (x&0x0001ffff) |
||||||
|
return -32767; |
||||||
|
else |
||||||
|
return 32767; |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
/*
|
||||||
|
K0 = 1 |
||||||
|
K1 = log(2) |
||||||
|
K2 = 3-4*log(2) |
||||||
|
K3 = 3*log(2) - 2 |
||||||
|
*/ |
||||||
|
#define D0 16384 |
||||||
|
#define D1 11356 |
||||||
|
#define D2 3726 |
||||||
|
#define D3 1301 |
||||||
|
/* Input in Q11 format, output in Q16 */ |
||||||
|
static inline spx_word32_t spx_exp2(spx_word16_t x) |
||||||
|
{ |
||||||
|
int integer; |
||||||
|
spx_word16_t frac; |
||||||
|
integer = SHR16(x,11); |
||||||
|
if (integer>14) |
||||||
|
return 0x7fffffff; |
||||||
|
else if (integer < -15) |
||||||
|
return 0; |
||||||
|
frac = SHL16(x-SHL16(integer,11),3); |
||||||
|
frac = ADD16(D0, MULT16_16_Q14(frac, ADD16(D1, MULT16_16_Q14(frac, ADD16(D2 , MULT16_16_Q14(D3,frac)))))); |
||||||
|
return VSHR32(EXTEND32(frac), -integer-2); |
||||||
|
} |
||||||
|
|
||||||
|
/* Input in Q11 format, output in Q16 */ |
||||||
|
static inline spx_word32_t spx_exp(spx_word16_t x) |
||||||
|
{ |
||||||
|
if (x>21290) |
||||||
|
return 0x7fffffff; |
||||||
|
else if (x<-21290) |
||||||
|
return 0; |
||||||
|
else |
||||||
|
return spx_exp2(MULT16_16_P14(23637,x)); |
||||||
|
} |
||||||
|
#define M1 32767 |
||||||
|
#define M2 -21 |
||||||
|
#define M3 -11943 |
||||||
|
#define M4 4936 |
||||||
|
|
||||||
|
static inline spx_word16_t spx_atan01(spx_word16_t x) |
||||||
|
{ |
||||||
|
return MULT16_16_P15(x, ADD32(M1, MULT16_16_P15(x, ADD32(M2, MULT16_16_P15(x, ADD32(M3, MULT16_16_P15(M4, x))))))); |
||||||
|
} |
||||||
|
|
||||||
|
#undef M1 |
||||||
|
#undef M2 |
||||||
|
#undef M3 |
||||||
|
#undef M4 |
||||||
|
|
||||||
|
/* Input in Q15, output in Q14 */ |
||||||
|
static inline spx_word16_t spx_atan(spx_word32_t x) |
||||||
|
{ |
||||||
|
if (x <= 32767) |
||||||
|
{ |
||||||
|
return SHR16(spx_atan01(x),1); |
||||||
|
} else { |
||||||
|
int e = spx_ilog2(x); |
||||||
|
if (e>=29) |
||||||
|
return 25736; |
||||||
|
x = DIV32_16(SHL32(EXTEND32(32767),29-e), EXTRACT16(SHR32(x, e-14))); |
||||||
|
return SUB16(25736, SHR16(spx_atan01(x),1)); |
||||||
|
} |
||||||
|
} |
||||||
|
#else |
||||||
|
|
||||||
|
#ifndef M_PI |
||||||
|
#define M_PI 3.14159265358979323846 /* pi */ |
||||||
|
#endif |
||||||
|
|
||||||
|
#define C1 0.9999932946f |
||||||
|
#define C2 -0.4999124376f |
||||||
|
#define C3 0.0414877472f |
||||||
|
#define C4 -0.0012712095f |
||||||
|
|
||||||
|
|
||||||
|
#define SPX_PI_2 1.5707963268 |
||||||
|
static inline spx_word16_t spx_cos(spx_word16_t x) |
||||||
|
{ |
||||||
|
if (x<SPX_PI_2) |
||||||
|
{ |
||||||
|
x *= x; |
||||||
|
return C1 + x*(C2+x*(C3+C4*x)); |
||||||
|
} else { |
||||||
|
x = M_PI-x; |
||||||
|
x *= x; |
||||||
|
return NEG16(C1 + x*(C2+x*(C3+C4*x))); |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
#endif |
||||||
@ -0,0 +1,169 @@ |
|||||||
|
/* Copyright (C) 2007 Jean-Marc Valin
|
||||||
|
|
||||||
|
File: os_support.h |
||||||
|
This is the (tiny) OS abstraction layer. Aside from math.h, this is the |
||||||
|
only place where system headers are allowed. |
||||||
|
|
||||||
|
Redistribution and use in source and binary forms, with or without |
||||||
|
modification, are permitted provided that the following conditions are |
||||||
|
met: |
||||||
|
|
||||||
|
1. Redistributions of source code must retain the above copyright notice, |
||||||
|
this list of conditions and the following disclaimer. |
||||||
|
|
||||||
|
2. Redistributions in binary form must reproduce the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer in the |
||||||
|
documentation and/or other materials provided with the distribution. |
||||||
|
|
||||||
|
3. The name of the author may not be used to endorse or promote products |
||||||
|
derived from this software without specific prior written permission. |
||||||
|
|
||||||
|
THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR |
||||||
|
IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES |
||||||
|
OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE |
||||||
|
DISCLAIMED. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, |
||||||
|
INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES |
||||||
|
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR |
||||||
|
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) |
||||||
|
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, |
||||||
|
STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN |
||||||
|
ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE |
||||||
|
POSSIBILITY OF SUCH DAMAGE. |
||||||
|
*/ |
||||||
|
|
||||||
|
#ifndef OS_SUPPORT_H |
||||||
|
#define OS_SUPPORT_H |
||||||
|
|
||||||
|
#include <string.h> |
||||||
|
#include <stdio.h> |
||||||
|
#include <stdlib.h> |
||||||
|
|
||||||
|
#ifdef HAVE_CONFIG_H |
||||||
|
#include "config.h" |
||||||
|
#endif |
||||||
|
#ifdef OS_SUPPORT_CUSTOM |
||||||
|
#include "os_support_custom.h" |
||||||
|
#endif |
||||||
|
|
||||||
|
/** Speex wrapper for calloc. To do your own dynamic allocation, all you need to do is replace this function, speex_realloc and speex_free
|
||||||
|
NOTE: speex_alloc needs to CLEAR THE MEMORY */ |
||||||
|
#ifndef OVERRIDE_SPEEX_ALLOC |
||||||
|
static inline void *speex_alloc (int size) |
||||||
|
{ |
||||||
|
/* WARNING: this is not equivalent to malloc(). If you want to use malloc()
|
||||||
|
or your own allocator, YOU NEED TO CLEAR THE MEMORY ALLOCATED. Otherwise |
||||||
|
you will experience strange bugs */ |
||||||
|
return calloc(size,1); |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
/** Same as speex_alloc, except that the area is only needed inside a Speex call (might cause problem with wideband though) */ |
||||||
|
#ifndef OVERRIDE_SPEEX_ALLOC_SCRATCH |
||||||
|
static inline void *speex_alloc_scratch (int size) |
||||||
|
{ |
||||||
|
/* Scratch space doesn't need to be cleared */ |
||||||
|
return calloc(size,1); |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
/** Speex wrapper for realloc. To do your own dynamic allocation, all you need to do is replace this function, speex_alloc and speex_free */ |
||||||
|
#ifndef OVERRIDE_SPEEX_REALLOC |
||||||
|
static inline void *speex_realloc (void *ptr, int size) |
||||||
|
{ |
||||||
|
return realloc(ptr, size); |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
/** Speex wrapper for calloc. To do your own dynamic allocation, all you need to do is replace this function, speex_realloc and speex_alloc */ |
||||||
|
#ifndef OVERRIDE_SPEEX_FREE |
||||||
|
static inline void speex_free (void *ptr) |
||||||
|
{ |
||||||
|
free(ptr); |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
/** Same as speex_free, except that the area is only needed inside a Speex call (might cause problem with wideband though) */ |
||||||
|
#ifndef OVERRIDE_SPEEX_FREE_SCRATCH |
||||||
|
static inline void speex_free_scratch (void *ptr) |
||||||
|
{ |
||||||
|
free(ptr); |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
/** Copy n elements from src to dst. The 0* term provides compile-time type checking */ |
||||||
|
#ifndef OVERRIDE_SPEEX_COPY |
||||||
|
#define SPEEX_COPY(dst, src, n) (memcpy((dst), (src), (n)*sizeof(*(dst)) + 0*((dst)-(src)) )) |
||||||
|
#endif |
||||||
|
|
||||||
|
/** Copy n elements from src to dst, allowing overlapping regions. The 0* term
|
||||||
|
provides compile-time type checking */ |
||||||
|
#ifndef OVERRIDE_SPEEX_MOVE |
||||||
|
#define SPEEX_MOVE(dst, src, n) (memmove((dst), (src), (n)*sizeof(*(dst)) + 0*((dst)-(src)) )) |
||||||
|
#endif |
||||||
|
|
||||||
|
/** For n elements worth of memory, set every byte to the value of c, starting at address dst */ |
||||||
|
#ifndef OVERRIDE_SPEEX_MEMSET |
||||||
|
#define SPEEX_MEMSET(dst, c, n) (memset((dst), (c), (n)*sizeof(*(dst)))) |
||||||
|
#endif |
||||||
|
|
||||||
|
|
||||||
|
#ifndef OVERRIDE_SPEEX_FATAL |
||||||
|
static inline void _speex_fatal(const char *str, const char *file, int line) |
||||||
|
{ |
||||||
|
fprintf (stderr, "Fatal (internal) error in %s, line %d: %s\n", file, line, str); |
||||||
|
exit(1); |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
#ifndef OVERRIDE_SPEEX_WARNING |
||||||
|
static inline void speex_warning(const char *str) |
||||||
|
{ |
||||||
|
#ifndef DISABLE_WARNINGS |
||||||
|
fprintf (stderr, "warning: %s\n", str); |
||||||
|
#endif |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
#ifndef OVERRIDE_SPEEX_WARNING_INT |
||||||
|
static inline void speex_warning_int(const char *str, int val) |
||||||
|
{ |
||||||
|
#ifndef DISABLE_WARNINGS |
||||||
|
fprintf (stderr, "warning: %s %d\n", str, val); |
||||||
|
#endif |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
#ifndef OVERRIDE_SPEEX_NOTIFY |
||||||
|
static inline void speex_notify(const char *str) |
||||||
|
{ |
||||||
|
#ifndef DISABLE_NOTIFICATIONS |
||||||
|
fprintf (stderr, "notification: %s\n", str); |
||||||
|
#endif |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
#ifndef OVERRIDE_SPEEX_PUTC |
||||||
|
/** Speex wrapper for putc */ |
||||||
|
static inline void _speex_putc(int ch, void *file) |
||||||
|
{ |
||||||
|
FILE *f = (FILE *)file; |
||||||
|
fprintf(f, "%c", ch); |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
#define speex_fatal(str) _speex_fatal(str, __FILE__, __LINE__); |
||||||
|
#define speex_assert(cond) {if (!(cond)) {speex_fatal("assertion failed: " #cond);}} |
||||||
|
|
||||||
|
#ifndef RELEASE |
||||||
|
static inline void print_vec(float *vec, int len, char *name) |
||||||
|
{ |
||||||
|
int i; |
||||||
|
printf ("%s ", name); |
||||||
|
for (i=0;i<len;i++) |
||||||
|
printf (" %f", vec[i]); |
||||||
|
printf ("\n"); |
||||||
|
} |
||||||
|
#endif |
||||||
|
|
||||||
|
#endif |
||||||
|
|
||||||
@ -0,0 +1,379 @@ |
|||||||
|
/* Copyright (C) 2005 Jean-Marc Valin */ |
||||||
|
/**
|
||||||
|
@file pseudofloat.h |
||||||
|
@brief Pseudo-floating point |
||||||
|
* This header file provides a lightweight floating point type for |
||||||
|
* use on fixed-point platforms when a large dynamic range is |
||||||
|
* required. The new type is not compatible with the 32-bit IEEE format, |
||||||
|
* it is not even remotely as accurate as 32-bit floats, and is not |
||||||
|
* even guaranteed to produce even remotely correct results for code |
||||||
|
* other than Speex. It makes all kinds of shortcuts that are acceptable |
||||||
|
* for Speex, but may not be acceptable for your application. You're |
||||||
|
* quite welcome to reuse this code and improve it, but don't assume |
||||||
|
* it works out of the box. Most likely, it doesn't. |
||||||
|
*/ |
||||||
|
/*
|
||||||
|
Redistribution and use in source and binary forms, with or without |
||||||
|
modification, are permitted provided that the following conditions |
||||||
|
are met: |
||||||
|
|
||||||
|
- Redistributions of source code must retain the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer. |
||||||
|
|
||||||
|
- Redistributions in binary form must reproduce the above copyright |
||||||
|
notice, this list of conditions and the following disclaimer in the |
||||||
|
documentation and/or other materials provided with the distribution. |
||||||
|
|
||||||
|
- Neither the name of the Xiph.org Foundation nor the names of its |
||||||
|
contributors may be used to endorse or promote products derived from |
||||||
|
this software without specific prior written permission. |
||||||
|
|
||||||
|
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS |
||||||
|
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT |
||||||
|
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR |
||||||
|
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE FOUNDATION OR |
||||||
|
CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, |
||||||
|
EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, |
||||||
|
PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR |
||||||
|
PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF |
||||||
|
LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING |
||||||
|
NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS |
||||||
|
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. |
||||||
|
*/ |
||||||
|
|
||||||
|
#ifndef PSEUDOFLOAT_H |
||||||
|
#define PSEUDOFLOAT_H |
||||||
|
|
||||||
|
#include "arch.h" |
||||||
|
#include "os_support.h" |
||||||
|
#include "math_approx.h" |
||||||
|
#include <math.h> |
||||||
|
|
||||||
|
#ifdef FIXED_POINT |
||||||
|
|
||||||
|
typedef struct { |
||||||
|
spx_int16_t m; |
||||||
|
spx_int16_t e; |
||||||
|
} spx_float_t; |
||||||
|
|
||||||
|
static const spx_float_t FLOAT_ZERO = {0,0}; |
||||||
|
static const spx_float_t FLOAT_ONE = {16384,-14}; |
||||||
|
static const spx_float_t FLOAT_HALF = {16384,-15}; |
||||||
|
|
||||||
|
#define MIN(a,b) ((a)<(b)?(a):(b)) |
||||||
|
static inline spx_float_t PSEUDOFLOAT(spx_int32_t x) |
||||||
|
{ |
||||||
|
int e=0; |
||||||
|
int sign=0; |
||||||
|
if (x<0) |
||||||
|
{ |
||||||
|
sign = 1; |
||||||
|
x = -x; |
||||||
|
} |
||||||
|
if (x==0) |
||||||
|
{ |
||||||
|
spx_float_t r = {0,0}; |
||||||
|
return r; |
||||||
|
} |
||||||
|
e = spx_ilog2(ABS32(x))-14; |
||||||
|
x = VSHR32(x, e); |
||||||
|
if (sign) |
||||||
|
{ |
||||||
|
spx_float_t r; |
||||||
|
r.m = -x; |
||||||
|
r.e = e; |
||||||
|
return r; |
||||||
|
} |
||||||
|
else |
||||||
|
{ |
||||||
|
spx_float_t r; |
||||||
|
r.m = x; |
||||||
|
r.e = e; |
||||||
|
return r; |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
|
||||||
|
static inline spx_float_t FLOAT_ADD(spx_float_t a, spx_float_t b) |
||||||
|
{ |
||||||
|
spx_float_t r; |
||||||
|
if (a.m==0) |
||||||
|
return b; |
||||||
|
else if (b.m==0) |
||||||
|
return a; |
||||||
|
if ((a).e > (b).e) |
||||||
|
{ |
||||||
|
r.m = ((a).m>>1) + ((b).m>>MIN(15,(a).e-(b).e+1)); |
||||||
|
r.e = (a).e+1; |
||||||
|
} |
||||||
|
else |
||||||
|
{ |
||||||
|
r.m = ((b).m>>1) + ((a).m>>MIN(15,(b).e-(a).e+1)); |
||||||
|
r.e = (b).e+1; |
||||||
|
} |
||||||
|
if (r.m>0) |
||||||
|
{ |
||||||
|
if (r.m<16384) |
||||||
|
{ |
||||||
|
r.m<<=1; |
||||||
|
r.e-=1; |
||||||
|
} |
||||||
|
} else { |
||||||
|
if (r.m>-16384) |
||||||
|
{ |
||||||
|
r.m<<=1; |
||||||
|
r.e-=1; |
||||||
|
} |
||||||
|
} |
||||||
|
/*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/ |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
static inline spx_float_t FLOAT_SUB(spx_float_t a, spx_float_t b) |
||||||
|
{ |
||||||
|
spx_float_t r; |
||||||
|
if (a.m==0) |
||||||
|
return b; |
||||||
|
else if (b.m==0) |
||||||
|
return a; |
||||||
|
if ((a).e > (b).e) |
||||||
|
{ |
||||||
|
r.m = ((a).m>>1) - ((b).m>>MIN(15,(a).e-(b).e+1)); |
||||||
|
r.e = (a).e+1; |
||||||
|
} |
||||||
|
else |
||||||
|
{ |
||||||
|
r.m = ((a).m>>MIN(15,(b).e-(a).e+1)) - ((b).m>>1); |
||||||
|
r.e = (b).e+1; |
||||||
|
} |
||||||
|
if (r.m>0) |
||||||
|
{ |
||||||
|
if (r.m<16384) |
||||||
|
{ |
||||||
|
r.m<<=1; |
||||||
|
r.e-=1; |
||||||
|
} |
||||||
|
} else { |
||||||
|
if (r.m>-16384) |
||||||
|
{ |
||||||
|
r.m<<=1; |
||||||
|
r.e-=1; |
||||||
|
} |
||||||
|
} |
||||||
|
/*printf ("%f + %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/ |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
static inline int FLOAT_LT(spx_float_t a, spx_float_t b) |
||||||
|
{ |
||||||
|
if (a.m==0) |
||||||
|
return b.m>0; |
||||||
|
else if (b.m==0) |
||||||
|
return a.m<0; |
||||||
|
if ((a).e > (b).e) |
||||||
|
return ((a).m>>1) < ((b).m>>MIN(15,(a).e-(b).e+1)); |
||||||
|
else |
||||||
|
return ((b).m>>1) > ((a).m>>MIN(15,(b).e-(a).e+1)); |
||||||
|
|
||||||
|
} |
||||||
|
|
||||||
|
static inline int FLOAT_GT(spx_float_t a, spx_float_t b) |
||||||
|
{ |
||||||
|
return FLOAT_LT(b,a); |
||||||
|
} |
||||||
|
|
||||||
|
static inline spx_float_t FLOAT_MULT(spx_float_t a, spx_float_t b) |
||||||
|
{ |
||||||
|
spx_float_t r; |
||||||
|
r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15); |
||||||
|
r.e = (a).e+(b).e+15; |
||||||
|
if (r.m>0) |
||||||
|
{ |
||||||
|
if (r.m<16384) |
||||||
|
{ |
||||||
|
r.m<<=1; |
||||||
|
r.e-=1; |
||||||
|
} |
||||||
|
} else { |
||||||
|
if (r.m>-16384) |
||||||
|
{ |
||||||
|
r.m<<=1; |
||||||
|
r.e-=1; |
||||||
|
} |
||||||
|
} |
||||||
|
/*printf ("%f * %f = %f\n", REALFLOAT(a), REALFLOAT(b), REALFLOAT(r));*/ |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
static inline spx_float_t FLOAT_AMULT(spx_float_t a, spx_float_t b) |
||||||
|
{ |
||||||
|
spx_float_t r; |
||||||
|
r.m = (spx_int16_t)((spx_int32_t)(a).m*(b).m>>15); |
||||||
|
r.e = (a).e+(b).e+15; |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
|
||||||
|
static inline spx_float_t FLOAT_SHL(spx_float_t a, int b) |
||||||
|
{ |
||||||
|
spx_float_t r; |
||||||
|
r.m = a.m; |
||||||
|
r.e = a.e+b; |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
static inline spx_int16_t FLOAT_EXTRACT16(spx_float_t a) |
||||||
|
{ |
||||||
|
if (a.e<0) |
||||||
|
return EXTRACT16((EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e); |
||||||
|
else |
||||||
|
return a.m<<a.e; |
||||||
|
} |
||||||
|
|
||||||
|
static inline spx_int32_t FLOAT_EXTRACT32(spx_float_t a) |
||||||
|
{ |
||||||
|
if (a.e<0) |
||||||
|
return (EXTEND32(a.m)+(EXTEND32(1)<<(-a.e-1)))>>-a.e; |
||||||
|
else |
||||||
|
return EXTEND32(a.m)<<a.e; |
||||||
|
} |
||||||
|
|
||||||
|
static inline spx_int32_t FLOAT_MUL32(spx_float_t a, spx_word32_t b) |
||||||
|
{ |
||||||
|
return VSHR32(MULT16_32_Q15(a.m, b),-a.e-15); |
||||||
|
} |
||||||
|
|
||||||
|
static inline spx_float_t FLOAT_MUL32U(spx_word32_t a, spx_word32_t b) |
||||||
|
{ |
||||||
|
int e1, e2; |
||||||
|
spx_float_t r; |
||||||
|
if (a==0 || b==0) |
||||||
|
{ |
||||||
|
return FLOAT_ZERO; |
||||||
|
} |
||||||
|
e1 = spx_ilog2(ABS32(a)); |
||||||
|
a = VSHR32(a, e1-14); |
||||||
|
e2 = spx_ilog2(ABS32(b)); |
||||||
|
b = VSHR32(b, e2-14); |
||||||
|
r.m = MULT16_16_Q15(a,b); |
||||||
|
r.e = e1+e2-13; |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
/* Do NOT attempt to divide by a negative number */ |
||||||
|
static inline spx_float_t FLOAT_DIV32_FLOAT(spx_word32_t a, spx_float_t b) |
||||||
|
{ |
||||||
|
int e=0; |
||||||
|
spx_float_t r; |
||||||
|
if (a==0) |
||||||
|
{ |
||||||
|
return FLOAT_ZERO; |
||||||
|
} |
||||||
|
e = spx_ilog2(ABS32(a))-spx_ilog2(b.m-1)-15; |
||||||
|
a = VSHR32(a, e); |
||||||
|
if (ABS32(a)>=SHL32(EXTEND32(b.m-1),15)) |
||||||
|
{ |
||||||
|
a >>= 1; |
||||||
|
e++; |
||||||
|
} |
||||||
|
r.m = DIV32_16(a,b.m); |
||||||
|
r.e = e-b.e; |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
|
||||||
|
/* Do NOT attempt to divide by a negative number */ |
||||||
|
static inline spx_float_t FLOAT_DIV32(spx_word32_t a, spx_word32_t b) |
||||||
|
{ |
||||||
|
int e0=0,e=0; |
||||||
|
spx_float_t r; |
||||||
|
if (a==0) |
||||||
|
{ |
||||||
|
return FLOAT_ZERO; |
||||||
|
} |
||||||
|
if (b>32767) |
||||||
|
{ |
||||||
|
e0 = spx_ilog2(b)-14; |
||||||
|
b = VSHR32(b, e0); |
||||||
|
e0 = -e0; |
||||||
|
} |
||||||
|
e = spx_ilog2(ABS32(a))-spx_ilog2(b-1)-15; |
||||||
|
a = VSHR32(a, e); |
||||||
|
if (ABS32(a)>=SHL32(EXTEND32(b-1),15)) |
||||||
|
{ |
||||||
|
a >>= 1; |
||||||
|
e++; |
||||||
|
} |
||||||
|
e += e0; |
||||||
|
r.m = DIV32_16(a,b); |
||||||
|
r.e = e; |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
/* Do NOT attempt to divide by a negative number */ |
||||||
|
static inline spx_float_t FLOAT_DIVU(spx_float_t a, spx_float_t b) |
||||||
|
{ |
||||||
|
int e=0; |
||||||
|
spx_int32_t num; |
||||||
|
spx_float_t r; |
||||||
|
if (b.m<=0) |
||||||
|
{ |
||||||
|
speex_warning_int("Attempted to divide by", b.m); |
||||||
|
return FLOAT_ONE; |
||||||
|
} |
||||||
|
num = a.m; |
||||||
|
a.m = ABS16(a.m); |
||||||
|
while (a.m >= b.m) |
||||||
|
{ |
||||||
|
e++; |
||||||
|
a.m >>= 1; |
||||||
|
} |
||||||
|
num = num << (15-e); |
||||||
|
r.m = DIV32_16(num,b.m); |
||||||
|
r.e = a.e-b.e-15+e; |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
static inline spx_float_t FLOAT_SQRT(spx_float_t a) |
||||||
|
{ |
||||||
|
spx_float_t r; |
||||||
|
spx_int32_t m; |
||||||
|
m = SHL32(EXTEND32(a.m), 14); |
||||||
|
r.e = a.e - 14; |
||||||
|
if (r.e & 1) |
||||||
|
{ |
||||||
|
r.e -= 1; |
||||||
|
m <<= 1; |
||||||
|
} |
||||||
|
r.e >>= 1; |
||||||
|
r.m = spx_sqrt(m); |
||||||
|
return r; |
||||||
|
} |
||||||
|
|
||||||
|
#else |
||||||
|
|
||||||
|
#define spx_float_t float |
||||||
|
#define FLOAT_ZERO 0.f |
||||||
|
#define FLOAT_ONE 1.f |
||||||
|
#define FLOAT_HALF 0.5f |
||||||
|
#define PSEUDOFLOAT(x) (x) |
||||||
|
#define FLOAT_MULT(a,b) ((a)*(b)) |
||||||
|
#define FLOAT_AMULT(a,b) ((a)*(b)) |
||||||
|
#define FLOAT_MUL32(a,b) ((a)*(b)) |
||||||
|
#define FLOAT_DIV32(a,b) ((a)/(b)) |
||||||
|
#define FLOAT_EXTRACT16(a) (a) |
||||||
|
#define FLOAT_EXTRACT32(a) (a) |
||||||
|
#define FLOAT_ADD(a,b) ((a)+(b)) |
||||||
|
#define FLOAT_SUB(a,b) ((a)-(b)) |
||||||
|
#define REALFLOAT(x) (x) |
||||||
|
#define FLOAT_DIV32_FLOAT(a,b) ((a)/(b)) |
||||||
|
#define FLOAT_MUL32U(a,b) ((a)*(b)) |
||||||
|
#define FLOAT_SHL(a,b) (a) |
||||||
|
#define FLOAT_LT(a,b) ((a)<(b)) |
||||||
|
#define FLOAT_GT(a,b) ((a)>(b)) |
||||||
|
#define FLOAT_DIVU(a,b) ((a)/(b)) |
||||||
|
#define FLOAT_SQRT(a) (spx_sqrt(a)) |
||||||
|
|
||||||
|
#endif |
||||||
|
|
||||||
|
#endif |
||||||
@ -0,0 +1,219 @@ |
|||||||
|
/* Тест обёртки эхоподавления lib/speex_aec (SpeexDSP mdf).
|
||||||
|
* |
||||||
|
* Проверяет: |
||||||
|
* - жизненный цикл и NULL/невалидные аргументы (create/destroy/reset/process); |
||||||
|
* - линию задержки: passthrough на старте, наполнение до delay_frames, учёт дрейфа; |
||||||
|
* - накопление захвата/рендера чанками произвольного размера (480/960); |
||||||
|
* - собственно подавление эха: синтетический far-end (шум) + near-end (тон), |
||||||
|
* после адаптации остаток эха должен быть заметно меньше исходного эха. |
||||||
|
*/ |
||||||
|
#include <stdio.h> |
||||||
|
#include <string.h> |
||||||
|
#include <math.h> |
||||||
|
|
||||||
|
#include "../lib/speex_aec.h" |
||||||
|
#include "../lib/debug_config.h" |
||||||
|
#include "../lib/mem.h" |
||||||
|
|
||||||
|
#define RATE 48000 |
||||||
|
#define FS 960 /* 20 мс */ |
||||||
|
#define FILTER 14400 /* хвост 300 мс */ |
||||||
|
#define DELAY_FRAMES 2 |
||||||
|
#define ECHO_GAIN 0.5 |
||||||
|
#define N_FRAMES 160 |
||||||
|
#define ADAPT_FRAMES 100 /* фаза адаптации: только far-end */ |
||||||
|
#define DT_START (ADAPT_FRAMES) |
||||||
|
#define DT_FRAMES 20 /* double-talk: речь + эхо */ |
||||||
|
#define MEAS_START (DT_START + DT_FRAMES) |
||||||
|
|
||||||
|
static int g_failures = 0; |
||||||
|
|
||||||
|
#define CHECK(cond, ...) do { \ |
||||||
|
if (!(cond)) { \
|
||||||
|
printf(" FAIL: " __VA_ARGS__); printf("\n"); \
|
||||||
|
g_failures++; \
|
||||||
|
} \
|
||||||
|
} while (0) |
||||||
|
|
||||||
|
/* Детерминированный «шум» far-end (LCG). */ |
||||||
|
static uint32_t g_rng = 0x12345678; |
||||||
|
static int16_t aec_noise(void) { |
||||||
|
g_rng = g_rng * 1664525u + 1013904223u; |
||||||
|
return (int16_t)(((g_rng >> 16) & 0x3fff) - 0x1fff); |
||||||
|
} |
||||||
|
|
||||||
|
static int test_lifecycle_and_null(void) { |
||||||
|
speex_aec_t* a; |
||||||
|
int16_t pcm[FS], out[FS]; |
||||||
|
|
||||||
|
printf("[speex_aec] lifecycle + NULL guards\n"); |
||||||
|
memset(pcm, 0, sizeof(pcm)); |
||||||
|
|
||||||
|
a = speex_aec_create(RATE, FS, FILTER, DELAY_FRAMES); |
||||||
|
CHECK(a != NULL, "create failed"); |
||||||
|
if (!a) return 1; |
||||||
|
|
||||||
|
CHECK(speex_aec_create(0, FS, FILTER, DELAY_FRAMES) == NULL, "create(rate=0) should fail"); |
||||||
|
CHECK(speex_aec_create(RATE, 0, FILTER, DELAY_FRAMES) == NULL, "create(frame=0) should fail"); |
||||||
|
CHECK(speex_aec_create(RATE, FS, FS - 1, DELAY_FRAMES) == NULL, "create(filter<frame) should fail"); |
||||||
|
CHECK(speex_aec_create(RATE, FS, FILTER, -3) != NULL, "create(delay=-3) should clamp, not fail"); |
||||||
|
|
||||||
|
CHECK(speex_aec_process_capture(NULL, pcm, FS, out) == 0, "process(NULL vad) should return 0"); |
||||||
|
CHECK(speex_aec_process_capture(a, NULL, FS, out) == 0, "process(NULL pcm) should return 0"); |
||||||
|
CHECK(speex_aec_process_capture(a, pcm, FS, NULL) == 0, "process(NULL out) should return 0"); |
||||||
|
CHECK(speex_aec_delay_fill(NULL) == 0, "delay_fill(NULL) should be 0"); |
||||||
|
|
||||||
|
speex_aec_feed_playback(NULL, pcm, FS); /* не должно падать */ |
||||||
|
speex_aec_reset(NULL); /* не должно падать */ |
||||||
|
speex_aec_destroy(NULL); /* не должно падать */ |
||||||
|
|
||||||
|
speex_aec_destroy(a); |
||||||
|
return 0; |
||||||
|
} |
||||||
|
|
||||||
|
static int test_delay_line(void) { |
||||||
|
speex_aec_t* a; |
||||||
|
int16_t in[FS], out[FS]; |
||||||
|
int16_t pcm[FS]; |
||||||
|
int i, n; |
||||||
|
|
||||||
|
printf("[speex_aec] delay line + chunk accumulation\n"); |
||||||
|
for (i = 0; i < FS; i++) pcm[i] = (int16_t)(i & 0x7ff); |
||||||
|
|
||||||
|
a = speex_aec_create(RATE, FS, FILTER, 3); |
||||||
|
CHECK(a != NULL, "create failed"); |
||||||
|
if (!a) return 1; |
||||||
|
|
||||||
|
CHECK(speex_aec_delay_fill(a) == 0, "initial fill=%d (expect 0)", speex_aec_delay_fill(a)); |
||||||
|
|
||||||
|
/* захват до рендера: passthrough (линия пуста) */ |
||||||
|
memset(in, 0x11, sizeof(in)); |
||||||
|
n = speex_aec_process_capture(a, in, FS, out); |
||||||
|
CHECK(n == FS, "process_capture returned %d (expect %d)", n, FS); |
||||||
|
CHECK(memcmp(in, out, sizeof(in)) == 0, "passthrough expected at empty line"); |
||||||
|
|
||||||
|
/* рендер чанками по 480 → накопление до кадра */ |
||||||
|
speex_aec_feed_playback(a, pcm, 480); |
||||||
|
CHECK(speex_aec_delay_fill(a) == 0, "fill=%d after 480 (expect 0)", speex_aec_delay_fill(a)); |
||||||
|
speex_aec_feed_playback(a, pcm + 480, 480); |
||||||
|
CHECK(speex_aec_delay_fill(a) == 1, "fill=%d after 960 (expect 1)", speex_aec_delay_fill(a)); |
||||||
|
|
||||||
|
/* захват чанками по 480 → 0, потом кадр */ |
||||||
|
n = speex_aec_process_capture(a, in, 480, out); |
||||||
|
CHECK(n == 0, "capture 480 should accumulate (returned %d)", n); |
||||||
|
n = speex_aec_process_capture(a, in + 480, 480, out); |
||||||
|
CHECK(n == FS, "capture 960 should emit frame (returned %d)", n); |
||||||
|
|
||||||
|
/* докормить рендер до глубины 3 */ |
||||||
|
speex_aec_feed_playback(a, pcm, FS); |
||||||
|
speex_aec_feed_playback(a, pcm, FS); |
||||||
|
CHECK(speex_aec_delay_fill(a) == 3, "fill=%d (expect 3)", speex_aec_delay_fill(a)); |
||||||
|
|
||||||
|
/* обработка кадра с полной линией уменьшает её на 1 */ |
||||||
|
n = speex_aec_process_capture(a, in, FS, out); |
||||||
|
CHECK(n == FS, "process_capture returned %d", n); |
||||||
|
CHECK(speex_aec_delay_fill(a) == 2, "fill=%d after capture (expect 2)", speex_aec_delay_fill(a)); |
||||||
|
|
||||||
|
speex_aec_reset(a); |
||||||
|
CHECK(speex_aec_delay_fill(a) == 0, "fill=%d after reset (expect 0)", speex_aec_delay_fill(a)); |
||||||
|
|
||||||
|
speex_aec_destroy(a); |
||||||
|
return 0; |
||||||
|
} |
||||||
|
|
||||||
|
static int test_echo_suppression(void) { |
||||||
|
speex_aec_t* a; |
||||||
|
int16_t* render_hist; |
||||||
|
int16_t speech[FS], cap[FS], out[FS]; |
||||||
|
int k, i, n; |
||||||
|
double echo_pow = 0.0, resid_pow = 0.0; |
||||||
|
double dt_speech_pow = 0.0, dt_out_pow = 0.0, dt_cross = 0.0; |
||||||
|
int measured = 0, dt_measured = 0; |
||||||
|
|
||||||
|
printf("[speex_aec] echo suppression (far-end noise, near-end tone)\n"); |
||||||
|
|
||||||
|
render_hist = (int16_t*)u_malloc((uint32_t)(N_FRAMES * FS) * sizeof(int16_t)); |
||||||
|
CHECK(render_hist != NULL, "OOM render_hist"); |
||||||
|
if (!render_hist) return 1; |
||||||
|
for (k = 0; k < N_FRAMES; k++) |
||||||
|
for (i = 0; i < FS; i++) render_hist[k * FS + i] = aec_noise(); |
||||||
|
|
||||||
|
a = speex_aec_create(RATE, FS, FILTER, DELAY_FRAMES); |
||||||
|
CHECK(a != NULL, "create failed"); |
||||||
|
if (!a) { u_free(render_hist); return 1; } |
||||||
|
|
||||||
|
for (k = 0; k < N_FRAMES; k++) { |
||||||
|
for (i = 0; i < FS; i++) { |
||||||
|
double t = (double)(k * FS + i) / RATE; |
||||||
|
speech[i] = (int16_t)(4000.0 * sin(2.0 * M_PI * 440.0 * t)); |
||||||
|
} |
||||||
|
|
||||||
|
/* фазы: [0,ADAPT) — только эхо (адаптация); [DT_START,MEAS_START) — речь+эхо;
|
||||||
|
* [MEAS_START,N) — только эхо (замер чистого подавления). */ |
||||||
|
int double_talk = (k >= DT_START && k < MEAS_START); |
||||||
|
for (i = 0; i < FS; i++) { |
||||||
|
int echo_idx = k - DELAY_FRAMES; |
||||||
|
int16_t echo = echo_idx >= 0 ? (int16_t)(ECHO_GAIN * render_hist[echo_idx * FS + i]) : 0; |
||||||
|
cap[i] = double_talk ? (int16_t)(speech[i] + echo) : echo; |
||||||
|
} |
||||||
|
|
||||||
|
n = speex_aec_process_capture(a, cap, FS, out); |
||||||
|
CHECK(n == FS, "process_capture returned %d (expect %d)", n, FS); |
||||||
|
speex_aec_feed_playback(a, render_hist + k * FS, FS); |
||||||
|
|
||||||
|
if (double_talk) { |
||||||
|
for (i = 0; i < FS; i++) { |
||||||
|
dt_speech_pow += (double)speech[i] * speech[i]; |
||||||
|
dt_out_pow += (double)out[i] * out[i]; |
||||||
|
dt_cross += (double)speech[i] * out[i]; |
||||||
|
} |
||||||
|
dt_measured += FS; |
||||||
|
} else if (k >= MEAS_START) { |
||||||
|
for (i = 0; i < FS; i++) { |
||||||
|
int echo_idx = k - DELAY_FRAMES; |
||||||
|
int16_t echo = echo_idx >= 0 ? (int16_t)(ECHO_GAIN * render_hist[echo_idx * FS + i]) : 0; |
||||||
|
echo_pow += (double)echo * echo; |
||||||
|
resid_pow += (double)out[i] * out[i]; |
||||||
|
} |
||||||
|
measured += FS; |
||||||
|
} |
||||||
|
} |
||||||
|
|
||||||
|
{ |
||||||
|
double echo_rms = sqrt(echo_pow / measured); |
||||||
|
double resid_rms = sqrt(resid_pow / measured); |
||||||
|
double db = echo_rms > 0 ? 20.0 * log10(echo_rms / (resid_rms > 1.0 ? resid_rms : 1.0)) : 0.0; |
||||||
|
printf(" echo-only suppression: echo_rms=%.1f resid_rms=%.1f => %.1f dB\n", echo_rms, resid_rms, db); |
||||||
|
CHECK(echo_rms > 100.0, "echo too weak to measure (echo_rms=%.1f)", echo_rms); |
||||||
|
CHECK(resid_pow < 0.05 * echo_pow, "echo not suppressed: resid=%.0f echo=%.0f", resid_pow, echo_pow); |
||||||
|
} |
||||||
|
{ |
||||||
|
double corr = sqrt(dt_speech_pow * dt_out_pow); |
||||||
|
double sim = corr > 0.0 ? dt_cross / corr : 0.0; |
||||||
|
printf(" double-talk speech preservation: correlation=%.3f\n", sim); |
||||||
|
CHECK(sim > 0.7, "near-end speech distorted: correlation=%.3f", sim); |
||||||
|
} |
||||||
|
|
||||||
|
speex_aec_destroy(a); |
||||||
|
u_free(render_hist); |
||||||
|
return 0; |
||||||
|
} |
||||||
|
|
||||||
|
int main(void) { |
||||||
|
debug_config_init(); |
||||||
|
debug_set_level(DEBUG_LEVEL_INFO); |
||||||
|
debug_set_category_level(DEBUG_CATEGORY_AEC, DEBUG_LEVEL_INFO); |
||||||
|
|
||||||
|
printf("SpeexDSP AEC test\n"); |
||||||
|
|
||||||
|
test_lifecycle_and_null(); |
||||||
|
test_delay_line(); |
||||||
|
test_echo_suppression(); |
||||||
|
|
||||||
|
if (g_failures == 0) { |
||||||
|
printf("TEST PASSED\n"); |
||||||
|
return 0; |
||||||
|
} |
||||||
|
printf("TEST FAILED (%d failures)\n", g_failures); |
||||||
|
return 1; |
||||||
|
} |
||||||
Loading…
Reference in new issue