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