Original FIR coeff script from LOFAR1 found by Andre Gunst. The...

Original FIR coeff script from LOFAR1 found by Andre Gunst. The pfs_coeff_final.m creates Coefficient_16KKaiser.dat, but these do differ slightly from git/hdl/libraries/dsp/filter/src/hex/Coeffs16384Kaiser-quant.dat that are used in LOFAR1.

applications/lofar2/model/pfs_coeff_final.m

+close all;
+clear all;
+
+fs =200*10^6;
+% fft_size;
+N=1024;
+% taps per frequency channel
+L=16;
+% passband ripple (in dB);
+r_pass=10^(0.5/20);
+% stopband ripple (in dB);
+r_stop=-89;
+% upsampling factor
+Q=16;
+% word size for coefficients
+W=16;
+
+% computed filter length
+M1=N*L/Q;
+% interpolated filter length
+M2=N*L;
+Nfft = M2;
+% compute filter
+c=fircls1(M1-1,Q/N,10^(r_pass/20),10^(r_stop/20),'trace');
+
+figure(30)
+Nfft =length(c);
+hold on
+spectrum_Fircls1= abs(fftshift(fft(c/max(c),Nfft)))/max(abs(fftshift(fft(c/max(c),Nfft))));
+plot([-Nfft/2+1:Nfft/2]/Nfft*fs,20*log10(spectrum_Fircls1),'r')
+hold off
+legend('Other candidate','filter given by fircls1')
+xlabel('f(Hz)')
+ylabel('Gain(dB)')
+title('Compirason of the prototype filter window before FFT interpolation');
+
+
+
+% use fourier interpolation to create final filter
+f1=fft(kaiser(length(c),1).'.*c);%hanning(length(c)).'.*
+f2=zeros(1,M2);
+% copy the lower frequency half. 
+n=0:M1/2;
+
+figure(31)
+subplot(2,1,1)
+plot(-length(f1)/2+1:length(f1)/2,20*log10(fftshift(abs(f1))))
+xlabel('f(Hz)')
+title('Gain of the filter after FFT')
+ylabel('Gain(dB)')
+subplot(2,1,2)
+plot(-length(f1)/2+1:length(f1)/2,unwrap(angle(fftshift(f1))))
+title('Phase of the filter after FFT')
+xlabel('f(Hz)')
+ylabel('Gain(dB)')
+
+f2(1+n)=f1(1+n);
+% to make the impulse response exactly symmetric, we need to do a phase
+% correction
+% to make the impulse response exactly symmetric, we need to do a phase correction
+f2(1+n)=f2(1+n).*exp(-sqrt(-1)*(Q-1)*pi*n*1/M2);
+% create the upper part of the spectrum from the lower part
+f2(M2-n)=conj(f2(2+n));
+
+figure(32)
+subplot(2,1,1)
+plot(-length(f2)/2+1:length(f2)/2,20*log10(fftshift(abs(f2))))
+xlabel('f(Hz)')
+title('Gain of the filter after FFT')
+ylabel('Gain(dB)')
+subplot(2,1,2)
+hold on
+plot(-length(f2)/2+1:length(f2)/2,unwrap(fftshift(angle(f2))))
+plot(-length(f2)/2+1:length(f2)/2,unwrap(fftshift(angle(f2))),'r')
+hold off
+title('Phase of the filter after FFT')
+xlabel('f(Hz)')
+ylabel('Gain(dB)')
+
+% back to time domain
+c2=real(ifft(f2));
+
+% quantize the coefficients
+c3=c2/max(c2);
+coeffs_quant = double(uencode(c3,16,max(c3),'signed'));
+% q=1/(2^(W-1));
+% c3=q*round(c3/q);
+% c3=c3*max(c2)/max(c3);
+
+fid = fopen('Coefficient_16KKaiser.dat','w');
+fprintf(fid,'%i\n', coeffs_quant);
+fclose(fid);
+fprintf('\n\nWARNING !\n')
+fprintf('This coefficient recordered in Coefficient_16K.dat are scaled to one. It is necessary then to multiply them with the number of bits needed for the quanization: 2^N_bits \n\n')
+
+
+% freqz(c3,1,M2*32,200*10^6);
+% q=1/(2^(W-1));
+% c3=q*round(c3/q);
+% c3=c3*max(c2)/max(c3);
+
+figure(2);
+n=0:M2-1;
+plot([c2',c2(M2-n)',c3']);
+title('impulse response');
+legend('unquantized','reversed','quantized');
+grid on;
+
+figure(2);
+hist(c2-c3,100);
+title('quantization error histogram');
+xlabel('quantization error');
+ylabel('relative probability');
+
+figure(10);
+freqz(c3,1,M2*4);
+title('frequency response');
+Nfft = 16384*8;
+c3 = c3/max(c3);
+%spectrum_fft = abs(fftshift(fft(c3/max(c3),Nfft)))/max(abs(fftshift(fft(c3/max(c3),Nfft))));
+spectrum = abs(fftshift(fft(c3/max(c3),Nfft)))/max(abs(fftshift(fft(c3/max(c3),Nfft))));
+frequency = [-Nfft/2+1:Nfft/2]*200*10^6/Nfft;
+plot(frequency(length(frequency)/2:end),20*log10(spectrum(length(frequency)/2:end)))
+grid on
+title('Gain of the filter window with quantized coefficients (16bits) using Fircls1 window with the FFT interpolation method: blackman window')
+xlabel('f(Hz)')
+ylabel('Gain(dB)')
+legend('FFT interpolation')