Filter analysis script for the polyphase low-pass pre-emphasis filter

doc/firpreemphasis.m

@@ -0,0 +1,81 @@
+
+
+% Filter analysis for combined LPF and preemphasis 
+% IIR reference material:  http://jontio.zapto.org/hda1/preempiir.pdf
+%
+% Tested on GNU Octave, should work on Matlab
+
+cutoff_freq=15700
+if 1
+  FIR_PHASES=32      % Values for PiFmRds
+  sample_rate=44100  
+  tau=75e-6          
+  delta=1.96e-6      
+else
+  FIR_PHASES=4       % Values used in reference material
+  sample_rate=48000  % 48000*4 = 192000
+  tau=50e-6          % to verify results against example values
+  delta=7.96e-6      
+end
+
+
+FIR_SIZE=1024
+
+taup=1/(2*sample_rate*FIR_PHASES)*cot(1/(2*tau*sample_rate*FIR_PHASES))
+deltap=1/(2*sample_rate*FIR_PHASES)*cot(1/(2*delta*sample_rate*FIR_PHASES))
+bp=sqrt(-(taup^2) + sqrt( (taup.^4) + 8*(taup^2)*(deltap^2) )  )/2
+ap=sqrt(2*(bp^2) + (taup^2) )
+a0=( 2*ap + 1/(sample_rate*FIR_PHASES) )/(2*bp + 1/(sample_rate*FIR_PHASES) )
+a1=(-2*ap + 1/(sample_rate*FIR_PHASES) )/(2*bp + 1/(sample_rate*FIR_PHASES) )
+b1=( 2*bp - 1/(sample_rate*FIR_PHASES) )/(2*bp + 1/(sample_rate*FIR_PHASES) )
+
+x=0;
+y=0;
+
+for j=1:FIR_SIZE
+  sincpos=j-((FIR_SIZE+1)/2);
+  if (sincpos==0)
+    firlowpass(j)=  2 * cutoff_freq / (sample_rate*FIR_PHASES) ;
+    disp('Oops: sincpos hit zero');
+  else
+    firlowpass(j)= sin(2 * pi * cutoff_freq * sincpos / (sample_rate*FIR_PHASES) ) / (pi * sincpos) ;
+  end
+  y=a0*firlowpass(j) + a1*x + b1*y;
+  x=firlowpass(j);
+  firpreemph(j)=y;
+  win(j)=(.54 - .46 * cos(2*pi * (j) / (FIR_SIZE))) * FIR_PHASES ; 
+
+end
+
+fullfilt=firpreemph.*win;
+
+hold off
+
+for j=1:FIR_PHASES
+  phasefreqs=linspace(0,sample_rate*FIR_PHASES*((FIR_SIZE/FIR_PHASES)-1)/FIR_SIZE,FIR_SIZE/FIR_PHASES);
+  phasedb=20*log10(abs(fft(fullfilt(j:FIR_PHASES:FIR_SIZE))));
+  % real data, so only show the positive frequencies
+  l1=plot( phasefreqs(1:(end/2)+1) , phasedb(1:(end/2+1)));
+  hold on
+
+  %plot( fullfilt(j:FIR_PHASES:FIR_SIZE) )
+end
+
+fulldb= 20*log10(abs(fft(fullfilt))/FIR_PHASES);
+fullfreqs=linspace(0,sample_rate*FIR_PHASES*(FIR_SIZE-1)/FIR_SIZE,FIR_SIZE);
+l2=plot( fullfreqs(1:(end/2+1)) ,fulldb(1:(end/2+1)),'r' );
+
+
+l3=plot( [19000 19000 ] , [0 -59] , 'g' ,'linewidth',2 ); % Pilot frequency 
+
+
+legend([l1,l2,l3],'Individual Phases','Full Filter Response','Pilot Tone')
+
+title('PiFmRds Polyphase Filter Response')
+xlim([0 60000]);
+xlabel('Hz');
+ylabel('dB');
+hold off
+
+% semilogx( linspace(0,sample_rate*FIR_PHASES*(FIR_SIZE-1)/FIR_SIZE,FIR_SIZE) ,20*log10(abs(fft(fullfilt))/FIR_PHASES) , 'r' )
+

src/fm_mpx.c

@@ -117,13 +117,13 @@ int fm_mpx_open(char *filename, size_t len) {
         }
     
         // Choose a cutoff frequency for the low-pass FIR filter
-        float cutoff_freq = 15000;
+        float cutoff_freq = 15700;
         if(in_samplerate/2 < cutoff_freq) cutoff_freq = in_samplerate/2 * .8;
    
 
         // Create the low-pass FIR filter, with pre-emphasis
         double window, firlowpass, firpreemph , sincpos;
-        double gain=FIR_PHASES/25.0;
+        double gain=FIR_PHASES/25.0; // Why??? Maybe gain adjustment for preemphais
 
         // IIR pre-emphasis filter
         // Reference material:    http://jontio.zapto.org/hda1/preempiir.pdf