--- /dev/null
+% 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))))
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])
//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.
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<DUC_M;j++,k++){
n = ((B1SMUL*n_1)>>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);
}
}
}
projects / freetel-svn-tracking.git / commitdiff