From: baobrien Date: Wed, 11 Mar 2015 20:22:50 +0000 (+0000) Subject: Work on 8Ksps->80Ksps upconversion and CIC filter X-Git-Url: http://git.whiteaudio.com/gitweb/?a=commitdiff_plain;h=a7e1800ad1828e38b372dc7214aa63e724353e29;p=freetel-svn-tracking.git Work on 8Ksps->80Ksps upconversion and CIC filter git-svn-id: https://svn.code.sf.net/p/freetel/code@2066 01035d8c-6547-0410-b346-abe4f91aad63 --- diff --git a/codec2-dev/octave/ciccomp.m b/codec2-dev/octave/ciccomp.m new file mode 100644 index 00000000..936a1013 --- /dev/null +++ b/codec2-dev/octave/ciccomp.m @@ -0,0 +1,70 @@ +% ciccomp.m +% Brady O'Brien 9 Mar 2015 +% CIC Filter compensation helper + +graphics_toolkit ("gnuplot"); + +cicn = 10; %delay for CIC filter +N = 10; %input interpolation rate +csf = 256; %scaling factor for CIC filter conversion +fd = 80e3; %DAC frequency +fs = fd/N; %Input sampling frequency +fi= fd; %freq of DSP +fc1 = fi/4; %Frequency of initial upconversion +ciccb=[-0.029626 0.252638 -2.304683 16.332166 -2.304683 0.252638 -0.029626]; +%ciccb = ciccb(1:2:length(ciccb)) +t = (1:fs); + +fdin = zeros(1,length(t)); +fdin(1)=1; +fdin=fdin+sin(pi*t*(2000/fs)); +%fdin = filter(b,1,fdin); + +figure(4) +plot(20*log10(abs(fft(fdin)))) +fdcout = zeros(1,length(t)); +fdin = int64(fdin*1024); + +combout1=0; +combout2=0; +combout3=0; +%combout4=0; +ccbuf1=zeros(1,cicn/N); +ccbuf2=zeros(1,cicn/N); +ccbuf3=zeros(1,cicn/N); +%ccbuf4=zeros(1,cicn); + +for ii=1:length(fdin) + combout1 = fdin(ii) - ccbuf1(end); + combout2 = combout1 - ccbuf2(end); + combout3 = combout2 - ccbuf3(end); + %combout4 = combout3 - ccbuf4(end); + ccbuf1(2:end)=ccbuf1(1:end-1); + ccbuf2(2:end)=ccbuf2(1:end-1); + ccbuf3(2:end)=ccbuf3(1:end-1); + %ccbuf4(2:end)=ccbuf4(1:end-1); + ccbuf1(1)=fdin(ii); + ccbuf2(1)=combout1; + ccbuf3(1)=combout2; + %ccbuf4(1)=combout3; + fdcout(ii)=combout3; +end + +intout1=0; +intout2=0; +intout3=0; +%intout4=0; +fdnext = zeros(1,length(t)*N); +fdnext(1:N:length(t)*N) = fdcout; %Interpolate +fdi1 = fdnext; + +for ii=1:length(fdnext) + intout1 = fdnext(ii) + intout1; + intout2 = intout1 + intout2; + intout3 = intout2 + intout3; + %intout4 = intout3 + intout4; + fdnext(ii)=intout3; +end + +figure(5) +plot(20*log10(abs(fft(fdnext)))) diff --git a/codec2-dev/octave/dacres.m b/codec2-dev/octave/dacres.m index 3d179639..34f7ee9d 100644 --- a/codec2-dev/octave/dacres.m +++ b/codec2-dev/octave/dacres.m @@ -5,52 +5,155 @@ graphics_toolkit ("gnuplot"); +cicn = 5; %delay for CIC filter +N = 5; %input interpolation rate +M = 25; %dac interpolation rate +csf = 1024; %scaling factor for CIC filter conversion +fd = 2E6; %DAC frequency +fs = 2E6/M/N; %Input sampling frequency +fi = 2E6/M; %freq of first interpolation +fb = 7E5; %Bandpass frequency +fc1 = fi/4; %Frequency of initial upconversion +ciccb=[-0.029626 0.252638 -2.304683 16.332166 -2.304683 0.252638 -0.029626]; %CIC Compensation FIR +pcicfb = fir1(20,.5); %Interpolation LPF Fir +s1fir = filter(ciccb,1,pcicfb); %Combine compensation and LPF + +t = (1:fs/2); +scpin = .1*sin(t*pi*2*(3000/8000))*i; % initial complex input +scpin(1) = 1+i; + +intstage1 = zeros(1,2*length(scpin)); %First stage of interpolation, 2x +intstage1(1:2:2*length(scpin))=scpin; +%scpin = filter(s1fir,1,scpin); + +scireal = int32(filter(s1fir,1,real(intstage1))*csf); %separate input into real and imiginary and apply pre-distortion +sciimag = int32(filter(s1fir,1,imag(intstage1))*csf); % also convert to integer. CIC integrator needs integer to work properly + + +%Apply 3 stage comb to real +fdin = scireal; +combout1=0; +combout2=0; +combout3=0; +combout4=0; +ccbuf1=zeros(1,cicn/N); +ccbuf2=zeros(1,cicn/N); +ccbuf3=zeros(1,cicn/N); +ccbuf4=zeros(1,cicn/N); + +for ii=1:length(fdin) + combout1 = fdin(ii) - ccbuf1(end); + combout2 = combout1 - ccbuf2(end); + combout3 = combout2 - ccbuf3(end); + combout4 = combout3 - ccbuf4(end); + ccbuf1(2:end)=ccbuf1(1:end-1); + ccbuf2(2:end)=ccbuf2(1:end-1); + ccbuf3(2:end)=ccbuf3(1:end-1); + ccbuf4(2:end)=ccbuf4(1:end-1); + ccbuf1(1)=fdin(ii); + ccbuf2(1)=combout1; + ccbuf3(1)=combout2; + ccbuf4(1)=combout3; + fdcout(ii)=combout4; +end + +intout1=0; +intout2=0; +intout3=0; +intout4=0; +fdnext = zeros(1,length(fdcout)*N); +fdnext(1:N:length(fdcout)*N) = fdcout; %Interpolate + +for ii=1:length(fdnext) + intout1 = fdnext(ii) + intout1; + intout2 = intout1 + intout2; + intout3 = intout2 + intout3; + intout4 = intout3 + intout4; + fdnext(ii)=intout4; +end +scoreal=single(fdnext/(2**20)); + +fdin=sciimag; +%Apply 3 stage comb to imag +combout1=0; +combout2=0; +combout3=0; +combout4=0; +ccbuf1=zeros(1,cicn/N); +ccbuf2=zeros(1,cicn/N); +ccbuf3=zeros(1,cicn/N); +ccbuf4=zeros(1,cicn/N); + +for ii=1:length(fdin) + combout1 = fdin(ii) - ccbuf1(end); + combout2 = combout1 - ccbuf2(end); + combout3 = combout2 - ccbuf3(end); + combout4 = combout3 - ccbuf4(end); + ccbuf1(2:end)=ccbuf1(1:end-1); + ccbuf2(2:end)=ccbuf2(1:end-1); + ccbuf3(2:end)=ccbuf3(1:end-1); + ccbuf4(2:end)=ccbuf4(1:end-1); + ccbuf1(1)=fdin(ii); + ccbuf2(1)=combout1; + ccbuf3(1)=combout2; + ccbuf4(1)=combout3; + fdcout(ii)=combout4; +end + +intout1=0; +intout2=0; +intout3=0; +intout4=0; +fdnext = zeros(1,length(fdcout)*N); +fdnext(1:N:length(fdcout)*N) = fdcout; %Interpolate + +for ii=1:length(fdnext) + intout1 = fdnext(ii) + intout1; + intout2 = intout1 + intout2; + intout3 = intout2 + intout3; + intout4 = intout3 + intout4; + fdnext(ii)=intout4; +end +scoimag=single(fdnext/(2**20)); + +%Convert complex to real and shift up to Fs/4 +ssout = scoreal.*cos(pi*.5*(1:length(scoreal))) + scoimag.*sin(pi*.5*(1:length(scoimag))); -M = 25; %interpolation -fs = 2E6/M; %Input sampling frequency -fb = 7E5; %Bandpass frequency -f1 = fs/4; -f2 = f1 + 4E3; -f3 = f1 - 4E3; t = (0:(fs-1)); beta1 = 0.999; -b1x = -2*sqrt(beta1)*cos(2*pi*fb/(fs*M)) -beta2 = beta1 - (1-0.999)*M; +b1x = -2*sqrt(beta1)*cos(2*pi*(fb/fd)) +beta2 = beta1 - (1-beta1)*M; -s1 = [fs zeros(1,fs-1)]; % noise floor, continuous interferers -s2 = 100*4*cos(t*2*pi*f2/fs); % wanted signal 40dB above interferers -s3 = 100*3*cos(t*2*pi*f3/fs); % another wanted frequency -s4 = 100*4*cos(t*2*pi*f1/fs); % wanted signal at center frequency -s = s1 + s2 + s3 + s4; +sducin = ssout; -s2 = filter([1 0 beta2],1,s); %pre interpolation notch filter to equalize bandpass after interpolation -s3 = zeros(1,length(s2)*M); %interpolate by zero-stuffing -s3(1:M:length(s2)*M) = s2; -s4 = filter(1,[1 b1x beta1],s3); %select wanted signal +sduceq = filter([1 0 beta2],1,sducin); %pre interpolation notch filter to equalize bandpass after interpolation +sducinterp = zeros(1,length(sduceq)*M); %interpolate by zero-stuffing +sducinterp(1:M:length(sduceq)*M) = sduceq; +sdac = filter(1,[1 b1x beta1],sducinterp); %select wanted signal figure(1) subplot(211) -plot(20*log10(abs(fft(s)/fs))) +plot(20*log10(abs(fft(sducin)/fi))) grid -axis([0 fs/2 -10 50]) +axis([0 fi/2 -20 50]) title('Output from modem'); subplot(212) -plot(20*log10(abs(fft(s2/fs)))) +plot(20*log10(abs(fft(sduceq/fi)))) grid -axis([0 fs/2 -10 70]) +axis([0 fi/2 -20 70]) title('After Pre-eq'); figure(2) subplot(211) -plot(20*log10(abs(fft(s3)/fs))) +plot(20*log10(abs(fft(sducinterp)/fd))) grid -axis([0 (fs/2)*M -10 80]) +axis([0 (fd/2) -20 80]) title('After interpolation'); subplot(212) -plot(20*log10(abs(fft(s4)/fs))) +plot(20*log10(abs(fft(sdac)/fd))) grid title('After bandpass'); -axis([0 (fs/2)*M -10 80]) +axis([0 (fd/2) -20 80]) diff --git a/codec2-dev/octave/octave-workspace b/codec2-dev/octave/octave-workspace index e69de29b..2f7f2224 100644 Binary files a/codec2-dev/octave/octave-workspace and b/codec2-dev/octave/octave-workspace differ diff --git a/codec2-dev/stm32/src/iir_duc.c b/codec2-dev/stm32/src/iir_duc.c index d4661131..4c69d0dc 100644 --- a/codec2-dev/stm32/src/iir_duc.c +++ b/codec2-dev/stm32/src/iir_duc.c @@ -47,7 +47,7 @@ //IIR and FIR filter states. Global for go fast. float f_1,f_2,f; -int n_1,n_2,n; +int n1_1,n1_2,n1,n2_1,n2_2,n2; /* Upconvert and bandpass filter a chunk of spectrum from Fs/M to Fs. We're going for 700khz here. @@ -64,17 +64,23 @@ void iir_upconv(float modin[], unsigned short dac_out[]){ f_2 = f_1; f_1 = modin[i]; //Scale fir output and convert to fixed. m = (int)((f/(IN_SCALE))*DAC_SCALE_2); //Scale fir output and convert to fixed - n = (m) + ((B1SMUL*n_1)>>B1SHFT) - ((B1MUL*n_2)>>B1SHFT); //Apply one cycle of IIR. This feeds the fir-ed sample into the output filter - n_2 = n_1; - n_1 = n; - dac_out[k]=(unsigned short)(n+DAC_SCALE_2); + n1 = m + ((B1SMUL*n1_1)>>B1SHFT) - ((B1MUL*n1_2)>>B1SHFT); //Apply one cycle of IIR. This feeds the fir-ed sample into the output filter + n1_2 = n1_1; + n1_1 = n1; + n2 = n1 + ((B1SMUL*n2_1)>>B1SHFT) - ((B1MUL*n2_2)>>B1SHFT); + n2_2 = n2_1; + n2_1 = n2; + dac_out[k]=(unsigned short)(n2+DAC_SCALE_2); k++; //now do the rest of the filtering. Because we're zero-stuffing we can neglect the sample from the fir filter. for(j=1;j>B1SHFT) - ((B1MUL*n_2)>>B1SHFT); n_2 = n_1; n_1 = n; - dac_out[k]=(unsigned short)((n)+DAC_SCALE_2); + n2 = n1 + ((B1SMUL*n2_1)>>B1SHFT) - ((B1MUL*n2_2)>>B1SHFT); + n2_2 = n2_1; + n2_1 = n2; + dac_out[k]=(unsigned short)((n2)+DAC_SCALE_2); } } }