I've actually seen this version doing as well as Fldigi in a simulation

grc_tests/bpsk31_tedvar.m

@@ -49,7 +49,7 @@ end
 
 function variance=run_var(snr, which_ted)
     disp('ran a command')
-    out_vect=shrun(sprintf('cat /tmp/psk31-raw-data | csdr awgn_cc %d --awgnfile /tmp/psk31-gaussian-noise | csdr timing_recovery_cc %s 256 0.5 2 --add_q --output_indexes | CSDR_FIXED_BUFSIZE=1048576 csdr normalized_timing_variance_u32_f 256 85', snr, which_ted), 'float32', 1);
+    out_vect=shrun(sprintf('cat /tmp/psk31-raw-data | csdr awgn_cc %d --awgnfile /tmp/psk31-gaussian-noise | csdr bandpass_fir_fft_cc $(csdr "=-31.25/8e3") $(csdr "=31.25/8e3") $(csdr "=31.25/8e3") | csdr simple_agc_cc 0.0001 0.5 | csdr timing_recovery_cc %s 256 0.5 2 --add_q --output_indexes | CSDR_FIXED_BUFSIZE=1048576 csdr normalized_timing_variance_u32_f 256 85', snr, which_ted), 'float32', 1);
     disp('run_var output:');
     out_vect'
     variance=out_vect(1);
@@ -77,7 +77,7 @@ function fmtplot(h)
     ylabel('Phase error variance [rad^2]');
 end
 
-snrs=-5:5:30
+snrs=-15:5:20
 %snrs=[10]
 error_values_gardner=mkvarplot('GARDNER',snrs);
 %{

grc_tests/test_bpsk_costas_loop_convertwavs.sh

@@ -1,9 +1,9 @@
 #!/bin/bash
 sox -r 48k -t f32 -c 2 /s/costas_nco -t wav -e floating-point /s/costas_nco.wav 
 sox -r 48k -t f32 -c 1 /s/costas_error -t wav -e floating-point /s/costas_error.wav 
-sox -r 48k -t f32 -c 1 /s/costas_dphase -t wav -e floating-point /s/costas_dphase.wav 
+sox -r 48k -t f32 -c 1 /s/costas_dphase -t wav -e floating-point --norm=-6 /s/costas_dphase.wav 
 sox -r 48k -t f32 -c 2 /s/costas_input -t wav -e floating-point /s/costas_input.wav 
 sox -r 48k -t f32 -c 2 /s/costas_output -t wav -e floating-point /s/costas_output.wav 
-
+ls -al /s/costas_nco.wav /s/costas_error.wav /s/costas_dphase.wav /s/costas_output.wav /s/costas_input.wav
 
 

libcsdr.c

@@ -2077,10 +2077,11 @@ void init_bpsk_costas_loop_cc(bpsk_costas_loop_state_t* s, int decision_directed
     //based on: http://gnuradio.squarespace.com/blog/2011/8/13/control-loop-gain-values.html
     float bandwidth_omega = 2*PI*bandwidth; //so that the bandwidth should be around 0.01 by default (2pi/100), and the damping_factor should be default 0.707
     float denomiator = 1+2*damping_factor*bandwidth_omega+bandwidth_omega*bandwidth_omega;
+    fprintf(stderr, "damp = %f, bw = %f, bwomega = %f\n", damping_factor, bandwidth, bandwidth_omega);
     s->alpha = (4*damping_factor*bandwidth_omega)/denomiator;
     s->beta = (4*bandwidth_omega*bandwidth_omega)/denomiator;
-    s->iir_temp = s->dphase = s->nco_phase = 0;
-    s->dphase_max=0.9*PI; //if it reached PI or -PI then it might actually hang and not come back
+    s->current_freq = s->dphase = s->nco_phase = 0;
+    s->dphase_max=bandwidth*PI; //this has been determined by experiment: if dphase is out of [-dphase_max; dphase_max] it might actually hang and not come back 
 }
 
 void bpsk_costas_loop_cc(complexf* input, complexf* output, int input_size, float* output_error, float* output_dphase, complexf* output_nco, bpsk_costas_loop_state_t* s)
@@ -2105,8 +2106,8 @@ void bpsk_costas_loop_cc(complexf* input, complexf* output, int input_size, floa
         }
         else error = PI*iof(output,i)*qof(output,i);
         if(output_error) output_error[i]=error;
-        s->dphase = error * s->alpha + s->iir_temp;
-        s->iir_temp += error * s->beta; //iir_temp could be named current_freq. See Tom Rondeau's article for better understanding.
+        s->current_freq += error * s->beta; 
+        s->dphase = error * s->alpha + s->current_freq;
         if(s->dphase>s->dphase_max) s->dphase=s->dphase_max;
         if(s->dphase<-s->dphase_max) s->dphase=-s->dphase_max;
         if(output_dphase) output_dphase[i]=s->dphase;
@@ -2187,7 +2188,7 @@ void simple_agc_cc(complexf* input, complexf* output, int input_size, float rate
         if(ideal_gain>max_gain) ideal_gain = max_gain;
         if(ideal_gain<=0) ideal_gain = 0;
         //*current_gain += (ideal_gain-(*current_gain))*rate;
-        *current_gain = (ideal_gain-(*current_gain))*rate + (*current_gain); //*rate_1minus;
+        *current_gain = (ideal_gain-(*current_gain))*rate + (*current_gain)*rate_1minus;
         //if(debugn<100) fprintf(stderr, "cgain: %g\n", *current_gain), debugn++;
         output[i].i=(*current_gain)*input[i].i;
         output[i].q=(*current_gain)*input[i].q;

libcsdr.h

@@ -367,7 +367,7 @@ typedef struct bpsk_costas_loop_state_s
     float alpha;
     float beta;
     int decision_directed;
-    float iir_temp;
+    float current_freq;
     float dphase;
     float nco_phase;
     float dphase_max;