csdr/libcsdr_gpl.c

263 lines
9.5 KiB
C
Raw Permalink Normal View History

2014-11-28 23:44:41 +08:00
/*
This file is part of libcsdr.
Copyright (c) Andras Retzler, HA7ILM <randras@sdr.hu>
Copyright (c) Warren Pratt, NR0V <warren@wpratt.com>
Copyright 2006,2010,2012 Free Software Foundation, Inc.
libcsdr is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
libcsdr is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with libcsdr. If not, see <http://www.gnu.org/licenses/>.
*/
#include "libcsdr_gpl.h"
#ifdef LIBCSDR_GPL
float shift_addition_cc(complexf *input, complexf* output, int input_size, shift_addition_data_t d, float starting_phase)
{
//The original idea was taken from wdsp:
//http://svn.tapr.org/repos_sdr_hpsdr/trunk/W5WC/PowerSDR_HPSDR_mRX_PS/Source/wdsp/shift.c
//However, this method introduces noise (from floating point rounding errors), which increases until the end of the buffer.
//fprintf(stderr, "cosd=%g sind=%g\n", d.cosdelta, d.sindelta);
float cosphi=cos(starting_phase);
float sinphi=sin(starting_phase);
float cosphi_last, sinphi_last;
for(int i=0;i<input_size;i++) //@shift_addition_cc: work
{
iof(output,i)=cosphi*iof(input,i)-sinphi*qof(input,i);
qof(output,i)=sinphi*iof(input,i)+cosphi*qof(input,i);
//using the trigonometric addition formulas
//cos(phi+delta)=cos(phi)cos(delta)-sin(phi)*sin(delta)
cosphi_last=cosphi;
sinphi_last=sinphi;
cosphi=cosphi_last*d.cosdelta-sinphi_last*d.sindelta;
sinphi=sinphi_last*d.cosdelta+cosphi_last*d.sindelta;
}
starting_phase+=d.rate*PI*input_size;
while(starting_phase>PI) starting_phase-=2*PI; //@shift_addition_cc: normalize starting_phase
while(starting_phase<-PI) starting_phase+=2*PI;
2014-11-28 23:44:41 +08:00
return starting_phase;
}
float shift_addition_fc(float *input, complexf* output, int input_size, shift_addition_data_t d, float starting_phase)
{
//The original idea was taken from wdsp:
//http://svn.tapr.org/repos_sdr_hpsdr/trunk/W5WC/PowerSDR_HPSDR_mRX_PS/Source/wdsp/shift.c
//However, this method introduces noise (from floating point rounding errors), which increases until the end of the buffer.
//fprintf(stderr, "cosd=%g sind=%g\n", d.cosdelta, d.sindelta);
float cosphi=cos(starting_phase);
float sinphi=sin(starting_phase);
float cosphi_last, sinphi_last;
for(int i=0;i<input_size;i++) //@shift_addition_cc: work
{
iof(output,i)=cosphi*input[i];
qof(output,i)=sinphi*input[i];
//using the trigonometric addition formulas
//cos(phi+delta)=cos(phi)cos(delta)-sin(phi)*sin(delta)
cosphi_last=cosphi;
sinphi_last=sinphi;
cosphi=cosphi_last*d.cosdelta-sinphi_last*d.sindelta;
sinphi=sinphi_last*d.cosdelta+cosphi_last*d.sindelta;
}
starting_phase+=d.rate*PI*input_size;
while(starting_phase>PI) starting_phase-=2*PI; //@shift_addition_cc: normalize starting_phase
while(starting_phase<-PI) starting_phase+=2*PI;
return starting_phase;
}
2014-11-28 23:44:41 +08:00
shift_addition_data_t shift_addition_init(float rate)
{
rate*=2;
shift_addition_data_t out;
out.sindelta=sin(rate*PI);
out.cosdelta=cos(rate*PI);
out.rate=rate;
return out;
}
#define SACCTEST_LOOPS 50
#define SACCTEST_STEP 10000
void shift_addition_cc_test(shift_addition_data_t d)
{
float phi=0;
float cosphi=cos(phi);
float sinphi=sin(phi);
float cosphi_last, sinphi_last;
int avg_size=(int)(2.0/d.rate+1.0); //average one period of sine
int avg_counter=0;
float avg=0;
printf("error_vector=[");
for(unsigned i=0;i<SACCTEST_STEP*SACCTEST_LOOPS;i++) //@shift_addition_cc: work
{
cosphi_last=cosphi;
sinphi_last=sinphi;
cosphi=cosphi_last*d.cosdelta-sinphi_last*d.sindelta;
sinphi=sinphi_last*d.cosdelta+cosphi_last*d.sindelta;
phi+=d.rate*PI;
while(phi>2*PI) phi-=2*PI; //@shift_addition_cc: normalize phase
if(i%SACCTEST_STEP==0)
2014-11-28 23:44:41 +08:00
{
avg_counter=avg_size;
avg=0;
}
if(avg_counter)
2014-11-28 23:44:41 +08:00
{
avg+=fabs(cosphi-cos(phi));
if(!--avg_counter) printf("%g ", avg/avg_size);
}
}
printf("]; error_vector_db=20*log10(error_vector); plot(error_vector_db);\n");
}
shift_addition_data_t decimating_shift_addition_init(float rate, int decimation)
{
return shift_addition_init(rate*decimation);
}
decimating_shift_addition_status_t decimating_shift_addition_cc(complexf *input, complexf* output, int input_size, shift_addition_data_t d, int decimation, decimating_shift_addition_status_t s)
{
//The original idea was taken from wdsp:
//http://svn.tapr.org/repos_sdr_hpsdr/trunk/W5WC/PowerSDR_HPSDR_mRX_PS/Source/wdsp/shift.c
//However, this method introduces noise (from floating point rounding errors), which increases until the end of the buffer.
//fprintf(stderr, "cosd=%g sind=%g\n", d.cosdelta, d.sindelta);
float cosphi=cos(s.starting_phase);
float sinphi=sin(s.starting_phase);
float cosphi_last, sinphi_last;
int i;
int k=0;
for(i=s.decimation_remain;i<input_size;i+=decimation) //@shift_addition_cc: work
{
iof(output,k)=cosphi*iof(input,i)-sinphi*qof(input,i);
qof(output,k)=sinphi*iof(input,i)+cosphi*qof(input,i);
k++;
//using the trigonometric addition formulas
//cos(phi+delta)=cos(phi)cos(delta)-sin(phi)*sin(delta)
cosphi_last=cosphi;
sinphi_last=sinphi;
cosphi=cosphi_last*d.cosdelta-sinphi_last*d.sindelta;
sinphi=sinphi_last*d.cosdelta+cosphi_last*d.sindelta;
}
s.decimation_remain=i-input_size;
s.starting_phase+=d.rate*PI*k;
s.output_size=k;
while(s.starting_phase>PI) s.starting_phase-=2*PI; //@shift_addition_cc: normalize starting_phase
while(s.starting_phase<-PI) s.starting_phase+=2*PI;
return s;
}
2014-11-28 23:44:41 +08:00
float agc_ff(float* input, float* output, int input_size, float reference, float attack_rate, float decay_rate, float max_gain, short hang_time, short attack_wait_time, float gain_filter_alpha, float last_gain)
{
/*
2014-12-03 19:02:21 +08:00
Notes on parameters (with some default values):
attack_rate = 0.01
decay_rate = 0.001
2014-11-28 23:44:41 +08:00
hang_time = (hang_time_ms / 1000) * sample_rate
hang_time is given in samples, and should be about 4ms.
hang_time can be switched off by setting it to zero (not recommended).
2014-11-28 23:44:41 +08:00
max_gain = pow(2, adc_bits)
max_gain should be no more than the dynamic range of your A/D converter.
gain_filter_alpha = 1 / ((fs/(2*PI*fc))+1)
2014-11-28 23:44:41 +08:00
>>> 1 / ((48000./(2*3.141592654*100))+1)
0.012920836043344543
>>> 1 / ((48000./(2*3.141592654*10))+1)
0.0013072857061786625
Literature:
2014-11-28 23:44:41 +08:00
ww.qsl.net/va3iul/Files/Automatic_Gain_Control.pdf
page 7 of http://www.arrl.org/files/file/Technology/tis/info/pdf/021112qex027.pdf
Examples:
http://svn.tapr.org/repos_sdr_hpsdr/trunk/W5WC/PowerSDR_HPSDR_mRX_PS/Source/wdsp/wcpAGC.c
GNU Radio's agc,agc2,agc3 have quite good ideas about this.
*/
register short hang_counter=0;
register short attack_wait_counter=0;
float gain=last_gain;
float last_peak=reference/last_gain; //approx.
float input_abs;
float error, dgain;
output[0]=last_gain*input[0]; //we skip this one sample, because it is easier this way
for(int i=1;i<input_size;i++) //@agc_ff
{
//The error is the difference between the required gain at the actual sample, and the previous gain value.
2014-11-28 23:44:41 +08:00
//We actually use an envelope detector.
input_abs=fabs(input[i]);
error=reference/input_abs-gain;
2014-11-28 23:44:41 +08:00
if(input[i]!=0) //We skip samples containing 0, as the gain would be infinity for those to keep up with the reference.
{
//An AGC is something nonlinear that's easier to implement in software:
//if the amplitude decreases, we increase the gain by minimizing the gain error by attack_rate.
//We also have a decay_rate that comes into consideration when the amplitude increases.
//The higher these rates are, the faster is the response of the AGC to amplitude changes.
//However, attack_rate should be higher than the decay_rate as we want to avoid clipping signals.
//that had a sudden increase in their amplitude.
//It's also important to note that this algorithm has an exponential gain ramp.
2014-11-28 23:44:41 +08:00
if(error<0) //INCREASE IN SIGNAL LEVEL
{
if(last_peak<input_abs)
2014-11-28 23:44:41 +08:00
{
2014-11-28 23:44:41 +08:00
attack_wait_counter=attack_wait_time;
last_peak=input_abs;
}
if(attack_wait_counter>0)
{
attack_wait_counter--;
//fprintf(stderr,"A");
2014-11-28 23:44:41 +08:00
dgain=0;
}
else
{
//If the signal level increases, we decrease the gain quite fast.
dgain=error*attack_rate;
2014-11-28 23:44:41 +08:00
//Before starting to increase the gain next time, we will be waiting until hang_time for sure.
hang_counter=hang_time;
2014-11-28 23:44:41 +08:00
}
}
else //DECREASE IN SIGNAL LEVEL
{
if(hang_counter>0) //Before starting to increase the gain, we will be waiting until hang_time.
{
hang_counter--;
dgain=0; //..until then, AGC is inactive and gain doesn't change.
}
else dgain=error*decay_rate; //If the signal level decreases, we increase the gain quite slowly.
}
2014-11-28 23:44:41 +08:00
gain=gain+dgain;
//fprintf(stderr,"g=%f dg=%f\n",gain,dgain);
}
if(gain>max_gain) gain=max_gain; //We also have to limit our gain, it can't be infinity.
if(gain<0) gain=0;
2014-11-28 23:44:41 +08:00
//output[i]=gain*input[i]; //Here we do the actual scaling of the samples.
//Here we do the actual scaling of the samples, but we run an IIR filter on the gain values:
output[i]=(gain=gain+last_gain-gain_filter_alpha*last_gain)*input[i]; //dc-pass-filter: freqz([1 -1],[1 -0.99]) y[i]=x[i]+y[i-1]-alpha*x[i-1]
2014-11-28 23:44:41 +08:00
//output[i]=input[i]*(last_gain+gain_filter_alpha*(gain-last_gain)); //LPF
last_gain=gain;
}
return gain; //this will be the last_gain next time
2014-11-28 23:44:41 +08:00
}
#endif