Subversion Repositories iqcorrection

% New IQ Correction algorithm, yet the third one

close all;

clear all;

e1 = 0.1; %gain error

%phase error

a1 = 10*pi/180;

%a1 = -10*pi/180;

%original noise

sgma = 0.01; %noise sigmage

%original sample length

%n_dat = 250000; % number of samples

n_dat = 100000; % number of samples

freq = 0.03;  % relative frequency

amplitude = 1;

%in phase signal, or I

x1=amplitude*cos(2*pi*(0:n_dat-1)*freq)+sgma*randn(1,n_dat);

%quadrature signal, or Q

y1=amplitude*(1+e1)*sin(2*pi*(0:n_dat-1)*freq+a1) + sgma*randn(1,n_dat);

reg_1 = 0;

reg_2 = 1;

%original mu

%mu_1 = 0.0002;

%mu_2 = 0.0001;

%mu adjusted to closest powers of two in order to match VHDL implementation

mu_1 = 0.000244;

mu_2 = 0.000122;

y2=zeros(1,n_dat);

y3=zeros(1,n_dat);

reg_1_sv = zeros(1,n_dat);

reg_2_sv = zeros(1,n_dat);

I_data = fopen('I_data_octave', 'w+t', 'native');

Q_data = fopen('Q_data_octave', 'w+t', 'native');

phase_estimate = fopen('phase_error_estimate_octave', 'w+t', 'native');

gain_estimate = fopen('gain_error_estimate_octave', 'w+t', 'native');

for nn=1:n_dat

    y2(nn) = y1(nn)-reg_1*x1(nn);

    reg_1_sv(nn) = reg_1;

    reg_1 = reg_1 + mu_1*x1(nn)*y2(nn);

    fprintf (I_data, '%f\r\n', x1(nn));

    fprintf (Q_data, '%f\r\n', y1(nn));

    fprintf (phase_estimate, '%f\r\n', reg_1);

    y3(nn) = y2(nn)*reg_2;

    reg_2_sv(nn) = reg_2;

    %original reg_2

    %reg_2 = reg_2+mu_2*(abs(x1(nn))^2 - abs(y3(nn))^2);

    %proposed reg_2

    reg_2 = reg_2+mu_2*((x1(nn))^2 - (y3(nn))^2);

    fprintf (gain_estimate, '%f\r\n', reg_2);

end

fclose(I_data);

fclose(Q_data);

fclose(phase_estimate);

fclose(gain_estimate);

figure(1)

subplot(2,1,1);

plot(reg_1_sv);

hold on;

plot([0 n_dat], [1 1]*a1*(1+e1), 'r');

hold off;

grid on;

title('Phase Error Estimate Trajectory')

subplot(2,1,2)

plot(reg_2_sv);

hold on

%original plot

%plot([0 n_dat], [1.0 1.0]/(1.0+e1), 'r');

%proposed plot

plot([0 n_dat], [1.0 1.0]/((1.0+e1)*(cos(a1))), 'r');

hold off

grid on

title('Gain Error Estimate Trajectory')

%added following line to make pngs

print("figure1_longer.png", "-dpng")

figure(2)

subplot(1,3,1)

plot(x1,y1)

hold on

plot(x1,y1,'r.')

plot([1 1 -1 -1 1],[1 -1 -1 1 1],'-r')

hold off

grid on

axis('square')

axis([-1.2 1.2 -1.2 1.2])

title('Lissajou Diagram, Gain & Phase Offset')

subplot(1,3,2)

plot(x1, y2)

hold on

plot(x1,y2, 'r.')

plot([1 1 -1 -1 1],[1 -1 -1 1 1],'-r')

hold off

grid on

axis('square')

axis([-1.2 1.2  -1.2 1.2])

title('Lissajou Diagram, Phase Offset Corrected')

subplot(1,3,3)

plot(x1, y3)

hold on

plot(x1,y3, 'r.')

plot([1 1 -1 -1 1],[1 -1 -1 1 1],'-r')

hold off

grid on

axis('square')

axis([-1.2 1.2  -1.2 1.2])

title('Lissajou Diagram, Gain & Phase Offset Corrected')

%added following line to make pngs

print("figure2_longer.png", "-dpng")

figure(3)

subplot(3,1,1)

%original MATLAB code used kaiser window function

%is commented out in this section.

%kaiser is in the package octave-windows

%this package turned out to be hard to install

%kaiser was replaced with blackman-harris window function.

%ww=kaiser(1024,12)';

ww=blackman(1024)';

ww=ww/sum(ww);

tt=n_dat-1023:n_dat;

plot(-0.5:1/1024:0.5-1/1024,fftshift(20*log10(abs(fft((x1(tt)+1i*y1(tt)).*ww)))))

% z=zeros(1,1024);

% for kk=1:1024

%    z(kk)=(x1(tt(kk))+1i*y1(tt(kk)))*ww(kk);

% end

% plot(-0.5:1/1024:0.5-1/1024,fftshift(20*log10(abs(fft(z)))))

grid on

axis([-0.5 0.5 -140 10])

title('Output Spectrum with Gain and Phase Offset')

subplot(3,1,2)

tt=n_dat-1023:n_dat;

% z=zeros(1,1024);

% for kk=1:1024

%     z(kk)=(x1(tt(kk))+1i*y2(tt(kk)))*ww(kk);

% end

% plot(-0.5:1/1024:0.5-1/1024,fftshift(20*log10(abs(fft(z)))));

plot(-0.5:1/1024:0.5-1/1024,fftshift(20*log10(abs(fft((x1(tt)+1i*y2(tt)).*ww)))))

grid on

axis([-0.5 0.5 -140 10])

title('Output Spectrum with Phase Offset Corrected')

subplot(3,1,3)

tt=n_dat-1023:n_dat;

% z=zeros(1,1024);

% for kk=1:1024

%     z(kk)=(x1(tt(kk))+1i*y3(tt(kk)))*ww(kk);

% end

% plot(-0.5:1/1024:0.5-1/1024,fftshift(20*log10(abs(fft(z)))));

plot(-0.5:1/1024:0.5-1/1024,fftshift(20*log10(abs(fft((x1(tt)+1i*y3(tt)).*ww)))))

grid on

axis([-0.5 0.5 -140 10])

title('Output Spectrum with Gain & Phase Offset Corrected')

%added following line to make pngs

print("figure3_longer.png", "-dpng")

disp('last reg_1 is ')

reg_1

disp('last reg_2 is ')

reg_2