00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039 #ifndef _VRAWGN_H_
00040 #define _VRAWGN_H_
00041
00042 #include <VrSigProc.h>
00043 #include <VrCycleCount.h>
00044
00045 #define IA 16807
00046 #define IM 2147483647
00047 #define AM (1.0/IM)
00048 #define IQ 127773
00049 #define IR 2836
00050 #define NTAB 32
00051 #define NDIV (1+(IM-1)/NTAB)
00052 #define EPS 1.2e-7
00053 #define RNMX (1.0-EPS)
00054
00055 template<class iType>
00056 class VrAWGN : public VrSigProc {
00057 protected:
00058 float maxAmp, psnr;
00059 float sigma;
00060 long seed;
00061 int iset;
00062 float gset;
00063 float ran1();
00064 float gasdev();
00065 public:
00066 virtual const char *name() { return "VrAWGN"; }
00067 virtual int work(VrSampleRange output, void *o[],
00068 VrSampleRange inputs[], void *i[]);
00069
00070
00071 VrAWGN(float snr = 72, float ma = 1024);
00072 void setSNR(float arg_snr);
00073 };
00074
00075
00076
00077
00078
00079 template<class iType> float
00080 VrAWGN<iType>::ran1()
00081 {
00082 int j;
00083 long k;
00084 static long iy=0;
00085 static long iv[NTAB];
00086 float temp;
00087
00088 if (seed <= 0 || !iy) {
00089 if (-seed < 1) seed=1;
00090 else seed = -seed;
00091 for (j=NTAB+7;j>=0;j--) {
00092 k=seed/IQ;
00093 seed=IA*(seed-k*IQ)-IR*k;
00094 if (seed < 0) seed += IM;
00095 if (j < NTAB) iv[j] = seed;
00096 }
00097 iy=iv[0];
00098 }
00099 k=(seed)/IQ;
00100 seed=IA*(seed-k*IQ)-IR*k;
00101 if (seed < 0)
00102 seed += IM;
00103 j=iy/NDIV;
00104 iy=iv[j];
00105 iv[j] = seed;
00106 temp=AM * iy;
00107 if (temp > RNMX)
00108 temp = RNMX;
00109 return temp;
00110 }
00111
00112
00113
00114
00115
00116 template<class iType> float
00117 VrAWGN<iType>::gasdev()
00118 {
00119 float fac,rsq,v1,v2;
00120
00121 iset = 1 - iset;
00122 if (iset) {
00123 do {
00124 v1=2.0*ran1()-1.0;
00125 v2=2.0*ran1()-1.0;
00126 rsq=v1*v1+v2*v2;
00127 } while (rsq >= 1.0 || rsq == 0.0);
00128 fac= sqrt(-2.0*log(rsq)/rsq) * sigma;
00129 gset=v1*fac;
00130 return v2*fac;
00131 }
00132 return gset;
00133 }
00134
00135 template<class iType> int
00136 VrAWGN<iType>::work (VrSampleRange output, void *ao[],
00137 VrSampleRange inputs[], void *ai[])
00138 {
00139 iType **i = (iType **)ai;
00140 iType **o = (iType **)ao;
00141 int size = output.size;
00142 while (size-- > 0) {
00143 iType temp = *i[0] + (iType)gasdev();
00144
00145 i[0]++;
00146 *o[0]++ = temp;
00147 }
00148 return output.size;
00149 }
00150
00151 template<class iType> void
00152 VrAWGN<iType>::setSNR(float arg_snr)
00153 {
00154 psnr = arg_snr;
00155 sigma = maxAmp / pow(10,(psnr/(double)20));
00156 printf ("snr %f ma %f sigma %f\n", psnr, maxAmp, sigma);
00157 }
00158
00159 template<class iType>
00160 VrAWGN<iType>::VrAWGN(float snr, float ma)
00161 : VrSigProc(1, sizeof(iType), sizeof(iType))
00162 {
00163 iset = 0;
00164 maxAmp = ma;
00165 setSNR (snr);
00166
00167
00168
00169 #if defined (__i386__)
00170
00171 seed = 3021;
00172
00173
00174
00175
00176
00177 #else
00178 seed = 3021;
00179 #endif
00180 }
00181 #endif