timidity_filter.cpp

Go to the documentation of this file.
00001 /*
00002 
00003     TiMidity -- Experimental MIDI to WAVE converter
00004     Copyright (C) 1995 Tuukka Toivonen <toivonen@clinet.fi>
00005 
00006     This program is free software; you can redistribute it and/or modify
00007     it under the terms of the GNU General Public License as published by
00008     the Free Software Foundation; either version 2 of the License, or
00009     (at your option) any later version.
00010 
00011     This program is distributed in the hope that it will be useful,
00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014     GNU General Public License for more details.
00015 
00016     You should have received a copy of the GNU General Public License
00017     along with this program; if not, write to the Free Software
00018     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
00019 
00020    filter.c: written by Vincent Pagel ( pagel@loria.fr ) 
00021  
00022    implements fir antialiasing filter : should help when setting sample
00023    rates as low as 8Khz.
00024 
00025    April 95
00026       - first draft
00027 
00028    22/5/95
00029       - modify "filter" so that it simulate leading and trailing 0 in the buffer
00030    */
00031 
00032 #include "pent_include.h"
00033 
00034 #ifdef USE_TIMIDITY_MIDI
00035 
00036 #include <cstdio>
00037 #include <cstring>
00038 #include <cmath>
00039 #include <cstdlib>
00040 #include "timidity.h"
00041 #include "timidity_common.h"
00042 #include "timidity_controls.h"
00043 #include "timidity_instrum.h"
00044 #include "timidity_filter.h"
00045 
00046 #ifdef NS_TIMIDITY
00047 namespace NS_TIMIDITY {
00048 #endif
00049 
00050 /*  bessel  function   */
00051 static float ino(float x)
00052 {
00053     float y, de, e, sde;
00054     int i;
00055     
00056     y = x / 2;
00057     e = 1.0;
00058     de = 1.0;
00059     i = 1;
00060     do {
00061         de = de * y / (float) i;
00062         sde = de * de;
00063         e += sde;
00064     } while (!( (e * 1.0e-08 - sde > 0) || (i++ > 25) ));
00065     return(e);
00066 }       
00067 
00068 /* Kaiser Window (symetric) */
00069 static void kaiser(float *w,int n,float beta)
00070 {
00071     float xind, xi;
00072     int i;
00073     
00074     xind = (float)((2*n - 1) * (2*n - 1));
00075     for (i =0; i<n ; i++) 
00076         {
00077             xi = (float)(i + 0.5);
00078             w[i] = ino((float)(beta * sqrt((double)(1. - 4 * xi * xi / xind))))
00079                 / ino((float)beta);
00080         }
00081 }
00082 
00083 /*
00084  * fir coef in g, cuttoff frequency in fc
00085  */
00086 static void designfir(float *g , float fc)
00087 {
00088     int i;
00089     float xi, omega, att, beta ;
00090     float w[ORDER2];
00091     
00092     for (i =0; i < ORDER2 ;i++) 
00093         {
00094             xi = (float) (i + 0.5);
00095             omega = (float)(PI * xi);
00096             g[i] = (float)(sin( (double) omega * fc) / omega);
00097         }
00098     
00099     att = 40.; /* attenuation  in  db */
00100     beta = (float) (exp(log((double)0.58417 * (att - 20.96)) * 0.4) + 0.07886 
00101         * (att - 20.96));
00102     kaiser( w, ORDER2, beta);
00103     
00104     /* Matrix product */
00105     for (i =0; i < ORDER2 ; i++)
00106         g[i] = g[i] * w[i];
00107 }
00108 
00109 /*
00110  * FIR filtering -> apply the filter given by coef[] to the data buffer
00111  * Note that we simulate leading and trailing 0 at the border of the 
00112  * data buffer
00113  */
00114 static void filter(sample_t *result,sample_t *data, sint32 length,float coef[])
00115 {
00116     sint32 sample,i,sample_window;
00117     sint16 peak = 0;
00118     float sum;
00119 
00120     /* Simulate leading 0 at the begining of the buffer */
00121      for (sample = 0; sample < ORDER2 ; sample++ )
00122         {
00123             sum = 0.0;
00124             sample_window= sample - ORDER2;
00125            
00126             for (i = 0; i < ORDER ;i++) 
00127                 sum += (float)(coef[i] *
00128                     ((sample_window<0)? 0.0 : data[sample_window++])) ;
00129             
00130             /* Saturation ??? */
00131             if (sum> 32767.) { sum=32767.; peak++; }
00132             if (sum< -32768.) { sum=-32768; peak++; }
00133             result[sample] = (sample_t) sum;
00134         }
00135 
00136     /* The core of the buffer  */
00137     for (sample = ORDER2; sample < length - ORDER + ORDER2 ; sample++ )
00138         {
00139             sum = 0.0;
00140             sample_window= sample - ORDER2;
00141 
00142             for (i = 0; i < ORDER ;i++) 
00143                 sum += data[sample_window++] * coef[i];
00144             
00145             /* Saturation ??? */
00146             if (sum> 32767.) { sum=32767.; peak++; }
00147             if (sum< -32768.) { sum=-32768; peak++; }
00148             result[sample] = (sample_t) sum;
00149         }
00150     
00151     /* Simulate 0 at the end of the buffer */
00152     for (sample = length - ORDER + ORDER2; sample < length ; sample++ )
00153         {
00154             sum = 0.0;
00155             sample_window= sample - ORDER2;
00156             
00157             for (i = 0; i < ORDER ;i++) 
00158                 sum += (float)(coef[i] *
00159                     ((sample_window>=length)? 0.0 : data[sample_window++])) ;
00160             
00161             /* Saturation ??? */
00162             if (sum> 32767.) { sum=32767.; peak++; }
00163             if (sum< -32768.) { sum=-32768; peak++; }
00164             result[sample] = (sample_t) sum;
00165         }
00166 
00167     if (peak)
00168         ctl->cmsg(CMSG_ERROR, VERB_NORMAL, 
00169                   "Saturation %2.3f %%.", 100.0*peak/ (float) length);
00170 }
00171 
00172 /***********************************************************************/
00173 /* Prevent aliasing by filtering any freq above the output_rate        */
00174 /*                                                                     */
00175 /* I don't worry about looping point -> they will remain soft if they  */
00176 /* were already                                                        */
00177 /***********************************************************************/
00178 void antialiasing(Sample *sp, sint32 output_rate )
00179 {
00180     sample_t *temp;
00181     int i;
00182     float fir_symetric[ORDER];
00183     float fir_coef[ORDER2];
00184     float freq_cut;  /* cutoff frequency [0..1.0] FREQ_CUT/SAMP_FREQ*/
00185  
00186 
00187     ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: Fsample=%iKHz",
00188               sp->sample_rate);
00189  
00190     /* No oversampling  */
00191     if (output_rate>=sp->sample_rate)
00192         return;
00193     
00194     freq_cut= (float) output_rate / (float) sp->sample_rate;
00195     ctl->cmsg(CMSG_INFO, VERB_NOISY, "Antialiasing: cutoff=%f%%",
00196               freq_cut*100.);
00197 
00198     designfir(fir_coef,freq_cut);
00199     
00200     /* Make the filter symetric */
00201     for (i = 0 ; i<ORDER2 ;i++) 
00202         fir_symetric[ORDER-1 - i] = fir_symetric[i] = fir_coef[ORDER2-1 - i];
00203     
00204     /* We apply the filter we have designed on a copy of the patch */
00205     temp = safe_Malloc<sample_t>(sp->data_length);
00206     memcpy(temp,sp->data,sp->data_length);
00207     
00208     filter(sp->data,temp,sp->data_length/sizeof(sample_t),fir_symetric);
00209     
00210     free(temp);
00211 }
00212 
00213 #ifdef NS_TIMIDITY
00214 };
00215 #endif
00216 
00217 #endif //USE_TIMIDITY_MIDI

Generated on Fri Jul 27 22:27:43 2007 for pentagram by  doxygen 1.4.7