/project/jni/sdl_mixer/timidity/filter.c

  1/*
  2
  3    TiMidity -- Experimental MIDI to WAVE converter
  4    Copyright (C) 1995 Tuukka Toivonen <toivonen@clinet.fi>
  5
  6    This program is free software; you can redistribute it and/or modify
  7    it under the terms of the GNU General Public License as published by
  8    the Free Software Foundation; either version 2 of the License, or
  9    (at your option) any later version.
 10
 11    This program is distributed in the hope that it will be useful,
 12    but WITHOUT ANY WARRANTY; without even the implied warranty of
 13    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 14    GNU General Public License for more details.
 15
 16    You should have received a copy of the GNU General Public License
 17    along with this program; if not, write to the Free Software
 18    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 19
 20   filter.c: written by Vincent Pagel ( pagel@loria.fr ) 
 21 
 22   implements fir antialiasing filter : should help when setting sample
 23   rates as low as 8Khz.
 24
 25   April 95
 26      - first draft
 27
 28   22/5/95
 29      - modify "filter" so that it simulate leading and trailing 0 in the buffer
 30   */
 31
 32#include <stdio.h>
 33#include <string.h>
 34#include <math.h>
 35#include <stdlib.h>
 36#include "config.h"
 37#include "common.h"
 38#include "ctrlmode.h"
 39#include "instrum.h"
 40#include "filter.h"
 41
 42/*  bessel  function   */
 43static float ino(float x)
 44{
 45    float y, de, e, sde;
 46    int i;
 47    
 48    y = x / 2;
 49    e = 1.0;
 50    de = 1.0;
 51    i = 1;
 52    do {
 53	de = de * y / (float) i;
 54	sde = de * de;
 55	e += sde;
 56    } while (!( (e * 1.0e-08 - sde > 0) || (i++ > 25) ));
 57    return(e);
 58}	
 59
 60/* Kaiser Window (symetric) */
 61static void kaiser(float *w,int n,float beta)
 62{
 63    float xind, xi;
 64    int i;
 65    
 66    xind = (float)((2*n - 1) * (2*n - 1));
 67    for (i =0; i<n ; i++) 
 68	{
 69	    xi = (float)(i + 0.5);
 70	    w[i] = ino((float)(beta * sqrt((double)(1. - 4 * xi * xi / xind))))
 71		/ ino((float)beta);
 72	}
 73}
 74
 75/*
 76 * fir coef in g, cuttoff frequency in fc
 77 */
 78static void designfir(float *g , float fc)
 79{
 80    int i;
 81    float xi, omega, att, beta ;
 82    float w[ORDER2];
 83    
 84    for (i =0; i < ORDER2 ;i++) 
 85	{
 86	    xi = (float) (i + 0.5);
 87	    omega = (float)(PI * xi);
 88	    g[i] = (float)(sin( (double) omega * fc) / omega);
 89	}
 90    
 91    att = 40.; /* attenuation  in  db */
 92    beta = (float) (exp(log((double)0.58417 * (att - 20.96)) * 0.4) + 0.07886 
 93	* (att - 20.96));
 94    kaiser( w, ORDER2, beta);
 95    
 96    /* Matrix product */
 97    for (i =0; i < ORDER2 ; i++)
 98	g[i] = g[i] * w[i];
 99}
100
101/*
102 * FIR filtering -> apply the filter given by coef[] to the data buffer
103 * Note that we simulate leading and trailing 0 at the border of the 
104 * data buffer
105 */
106static void filter(sample_t *result,sample_t *data, int32 length,float coef[])
107{
108    int32 sample,i,sample_window;
109    int16 peak = 0;
110    float sum;
111
112    /* Simulate leading 0 at the begining of the buffer */
113     for (sample = 0; sample < ORDER2 ; sample++ )
114	{
115	    sum = 0.0;
116	    sample_window= sample - ORDER2;
117	   
118	    for (i = 0; i < ORDER ;i++) 
119		sum += (float)(coef[i] *
120		    ((sample_window<0)? 0.0 : data[sample_window++])) ;
121	    
122	    /* Saturation ??? */
123	    if (sum> 32767.) { sum=32767.; peak++; }
124	    if (sum< -32768.) { sum=-32768; peak++; }
125	    result[sample] = (sample_t) sum;
126	}
127
128    /* The core of the buffer  */
129    for (sample = ORDER2; sample < length - ORDER + ORDER2 ; sample++ )
130	{
131	    sum = 0.0;
132	    sample_window= sample - ORDER2;
133
134	    for (i = 0; i < ORDER ;i++) 
135		sum += data[sample_window++] * coef[i];
136	    
137	    /* Saturation ??? */
138	    if (sum> 32767.) { sum=32767.; peak++; }
139	    if (sum< -32768.) { sum=-32768; peak++; }
140	    result[sample] = (sample_t) sum;
141	}
142    
143    /* Simulate 0 at the end of the buffer */
144    for (sample = length - ORDER + ORDER2; sample < length ; sample++ )
145	{
146	    sum = 0.0;
147	    sample_window= sample - ORDER2;
148	    
149	    for (i = 0; i < ORDER ;i++) 
150		sum += (float)(coef[i] *
151		    ((sample_window>=length)? 0.0 : data[sample_window++])) ;
152	    
153	    /* Saturation ??? */
154	    if (sum> 32767.) { sum=32767.; peak++; }
155	    if (sum< -32768.) { sum=-32768; peak++; }
156	    result[sample] = (sample_t) sum;
157	}
158
159    if (peak)
160	ctl->cmsg(CMSG_ERROR, VERB_NORMAL, 
161		  "Saturation %2.3f %%.", 100.0*peak/ (float) length);
162}
163
164/***********************************************************************/
165/* Prevent aliasing by filtering any freq above the output_rate        */
166/*                                                                     */
167/* I don't worry about looping point -> they will remain soft if they  */
168/* were already                                                        */
169/***********************************************************************/
170void antialiasing(Sample *sp, int32 output_rate )
171{
172    sample_t *temp;
173    int i;
174    float fir_symetric[ORDER];
175    float fir_coef[ORDER2];
176    float freq_cut;  /* cutoff frequency [0..1.0] FREQ_CUT/SAMP_FREQ*/
177 
178
179    ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: Fsample=%iKHz",
180	      sp->sample_rate);
181 
182    /* No oversampling  */
183    if (output_rate>=sp->sample_rate)
184	return;
185    
186    freq_cut= (float) output_rate / (float) sp->sample_rate;
187    ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: cutoff=%f%%",
188	      freq_cut*100.);
189
190    designfir(fir_coef,freq_cut);
191    
192    /* Make the filter symetric */
193    for (i = 0 ; i<ORDER2 ;i++) 
194	fir_symetric[ORDER-1 - i] = fir_symetric[i] = fir_coef[ORDER2-1 - i];
195    
196    /* We apply the filter we have designed on a copy of the patch */
197    temp = safe_malloc(sp->data_length);
198    memcpy(temp,sp->data,sp->data_length);
199    
200    filter(sp->data,temp,sp->data_length/sizeof(sample_t),fir_symetric);
201    
202    free(temp);
203}