Circular Convolution

Kadhiem Ayob March 3, 2012 Coded in Matlab
%circular convolution
%for testing you may use:
h = fir1(20,.3);
x = randn(1,1024);

%function y = conv_circ(h,x)

y = conv(h,x);

L1 = length(h);
L2 = length(x);

%add end to start, add start to end
temp = y(1:L1-1);
y(1:L1-1) = y(1:L1-1) + y(L2+(1:L1-1));
y(L2+(1:L1-1)) = y(L2+(1:L1-1)) + temp;

%compare to direct convolution
y2 = conv(h,x);
plot(y,'o-');hold
plot(y2,'r.--')
legend('circular','direct')

Power Computation of a Digital Stream

Kadhiem Ayob February 26, 2012 Coded in Matlab
clear all;

%example vector having ac & dc power
x = complex(randn(1,2048)+.2,randn(1,2048)-.1);

%total power, three equivalent methods
pwr1_total = mean(abs(x).^2);                %mean of squared values
pwr2_total = mean(real(x).^2 + imag(x).^2);
pwr3_total = mean(x .* conj(x));

%total expressed as rms
rms_total = sqrt(pwr1_total);

%dc power
pwr1_dc = mean(real(x))^2 + mean(imag(x))^2; %square of mean of values

%ac power
pwr1_ac = pwr1_total - pwr1_dc;              %mean of squares - square of mean

%relation of ac power to statistical variance 
pwr2_ac = var(x);              %approximately

%ac expressed as rms
rms1_ac = sqrt(pwr1_ac);

%ac relation to standard of deviation, std = sqrt(var)
rms2_ac = std(x);              %approximately

%dc relation to variance
pwr2_dc = pwr1_total - var(x); %approximately

fprintf('----------------------------------------------------\r');
fprintf('power(total),          : %1.5f, %1.5f, %1.5f\r',pwr1_total,pwr2_total,pwr3_total);
fprintf('rms(total)             : %1.5f\r',rms_total);
fprintf('power(ac),  variance   : %1.5f, %1.5f\r',pwr1_ac,pwr2_ac);
fprintf('rms(ac),    std        : %1.5f, %1.5f\r',rms1_ac,rms2_ac);
fprintf('power(dc), (total-var) : %1.5f, %1.5f\r',pwr1_dc,pwr2_dc);

Reading text files in Matlab

Kadhiem Ayob February 13, 2012 Coded in Matlab
clear all;

% create test file
data = [1 5 2 4 3 3 4 2 5 1];
filename = 'test_file.txt';
fid = fopen(filename, 'w');
fprintf(fid, '%d %d\n', data);
fclose(fid);

% (1) read file using fscanf 
fid = fopen(filename, 'r'); 
y1 = fscanf(fid, '%d\n'); %interleaves columns
fclose(fid);

% (2) read file using textread (or textscan)
[ya,yb] = textread(filename,'%d%d');
y2 = [ya yb];

% (3) read file using importdata 
y3 = importdata(filename);

% (4) read file using load 
y4 = load(filename);

disp('-------------------------------------------------------------')
disp('    original vector data')
disp(data)
disp('    file content using fprintf')
disp(y2)
disp('    vector created by fscanf')
disp(y1)
disp('     matrix created by:')
disp('     textread   importdata    load')
disp([y2 y3 y4])

Ideal interpolation filter

Kadhiem Ayob January 21, 2012 Coded in Matlab
clear all; close all;

%example bandlimited random input & parameters
x = filter(fir1(70,.1),1,randn(1,1024));

up = 3;             %Interpolation factor
cutoff = .3;

intorder = 6;
h1 = intfilt(up, intorder, cutoff); %ideal filter
h1 = up*h1/sum(h1);

h2 = fir1(2*up*intorder-2,cutoff);  %ordinary LPF 
h2 = up*h2/sum(h2);

%upsample
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;

x_f1 = filter(h1,1,x_up);
x_f2 = filter(h2,1,x_up);

figure;
subplot(3,1,1);hold
plot(x_f1,'o--');
plot(x_f2,'r.-');
legend('ideal output','fir1 output');

subplot(3,1,2);
plot(x_f1-x_f2);
legend('error');

subplot(3,1,3);hold
plot(h1,'.-');
plot(h2,'r.-');
legend('ideal filter','fir1 filter');

checking resampling in time domain

Kadhiem Ayob January 21, 2012 Coded in Matlab
clear all; close all;

%example bandlimited random input & parameters
x = filter(fir1(70,.1),1,randn(1,512));
h = fir1(30,.3);    %filter used for upsampling
up = 3;             %Interpolation factor
dn = 2;             %Decimation factor

%%%%%%%%%%%%%%%%%%%%% up/dn model %%%%%%%%%%%%%%%%%%%%%%
%upsample input
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;
x_f = filter(up*h,1,x_up);

%downsample signal by decimation
x_dn = x_f(1:dn:end);

delay = 30/2+1;  
figure;
subplot(2,1,1);hold;
plot(x_up,'o--');
plot(x_f(delay:end),'r.-');
legend('input with zeros','upsampled stage');

subplot(2,1,2);hold;
plot(x_up(1:dn:end),'o--');
plot(x_dn(ceil(delay/dn):end),'r.-');
legend('input with zeros','final signal');

Resampling filter performance

Kadhiem Ayob January 14, 20123 comments Coded in Matlab
%%%%%%%%%%%%%%%%%% inputs for model %%%%%%%%%%%%%%%%
clear all; close all;

%example bandlimited impulse input & parameters
x = filter(fir1(70,.1),1,[1 zeros(1,2^15-1)]);
Fs = 120;           %MHz original sample rate             
h = fir1(30,.3);    %filter used for upsampling
up = 3;             %Interpolation factor
dn = 2;             %Decimation factor
Fc = 12;             %MHz band centre (-Fs/2 ~ +Fs/2)
Fch = 0;             %MHz filter centre (-Fs*up/2 ~ +Fs*up/2)

%move signal to its centre
x = x.*exp(j*2*pi*(0:length(x)-1)*Fc/Fs);

%shift filter 
h = h.*exp(j*2*pi*(0:length(h)-1)*Fch/(Fs*up));

%%%%%%%%%%%%%%%%%%%%% model %%%%%%%%%%%%%%%%%%%%%%
%check signal in upsampled domain
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;
[P, F] = pwelch(x_up, [], 0, 2^16, Fs*up,'twosided'); 
F = F - max(F)/2; 
P = fftshift(P);
y(find(P == 0)) = -100;  %avoid log of zero 
y(find(P ~= 0)) = 10*log10(P(find(P ~= 0)));
P_dB = y - 10*log10(max(P));  %normalise

%check filter response in upsampled domain
H = fftshift(20*log10(abs(fft(h,2^16))));

subplot(2,1,1);
hold;grid;
plot(F, P_dB,'.-'); 
plot(F,H,'m--');
axis([min(F)-1 max(F)+1 -80 1]);
legend('upsampled signal','upsampling filter');

%check signal in downsampled domain
x_f = filter(h,1,x_up);
x_dn = x_f(1:dn:end);
[P, F] = pwelch(x_dn, [], 0, 2^16, Fs*up/dn,'twosided'); 
F = F - max(F)/2;
P = fftshift(P);
y(find(P == 0)) = -100;   %avoid log of zero 
y(find(P ~= 0)) = 10*log10(P(find(P ~= 0)));
P_dB = y - 10*log10(max(P)); %normalise

subplot(2,1,2)
plot(F,P_dB,'r.-');
grid;
axis([min(F)-1 max(F)+1 -80 1]);
legend('downsampled signal');

16QAM Modem model (Basic)

Kadhiem Ayob December 12, 2011 Coded in Matlab
clear all; close all;

n = 2^16;   %number of symbols
Fsym = 12.5; %Msps
Fs = 100;   %MHz, IF sampling frequency
Fc = 20;    %MHz, upconverter frequency

%generate 16QAM raw symbols
alphabet = [-1 -1/3 +1/3 +1];
tx = complex(randsrc(1,n,alphabet),randsrc(1,n,alphabet));

%pulse shaping, root raised cosine
h = firrcos(50,Fsym/4,.15,Fsym,'rolloff','sqrt');
tx_up1 = zeros(1,2*n);
tx_up1(1:2:end-1) = tx;
tx_shaped = filter(2*h,1,tx_up1);

%further upsampling by 4
tx_up2 = resample(tx_shaped,4,1);

%upconvert to 10MHz centre
f1 = exp(1i*2*pi*(0:8*n-1)*Fc/Fs);
tx_upconverted = tx_up2 .* f1;

%channel signal
rx_real = real(tx_upconverted);

%Rx shifts signal back to zero 
f2 = exp(1i*2*pi*(0:8*n-1)*-Fc/Fs);
rx_downconverted = rx_real.*f2;

%Rx downsamples back by 4
rx_dn1 = resample(rx_downconverted,1,4);

%rx matched filter
rx_dn2 = filter(2*h,1,rx_dn1);
rx_dn2 = rx_dn2(50:end); %remove initial zeros

%one phase should be correct(odd/even)
rx_e = rx_dn2(1:2:end); %odd
rx_o = rx_dn2(2:2:end); %even

subplot(3,1,1);plot(real(tx),imag(tx),'o');
axis([-2 +2 -2 +2]);grid;
title('Tx constellations')
subplot(3,1,2);plot(real(rx_e),imag(rx_e),'o');
axis([-2 +2 -2 +2]);grid;
title('Rx constellations, first phase')
subplot(3,1,3);plot(real(rx_o),imag(rx_o),'o');
axis([-2 +2 -2 +2]);grid;
title('Rx constellations, second phase')

Phase drift of NCO

Kadhiem Ayob December 10, 2011 Coded in Matlab
clear all; close all;

n = 50000;                            %number of test samples
Fs = 100;                             %sampling rate in Msps
Fo = 1.17;                            %target NCO frequency in MHz
M = [2^16 10000];                     %NCO phase resolution, two cases
for i = 1:2
    inc = round(M(i)*Fo/Fs);                         %phase increment
    lut = round(2047*cos(2*pi*(0:M(i)-1)/M(i)));     %lut, one cycle cosine data
    addr = 0:inc:inc*(n-1);
    addr = round(mod(addr,M(i)));
    addr(addr >= M(i)) = addr(addr >= M(i)) - M(i);  %check address overflow
    y(i,1:n) = lut(addr+1);                          %add 1 for matlab LUT
end

fprintf('increment value = %2.6f\r',M*Fo/Fs);

subplot(2,1,1);hold;
plot([y(1,1:1000) zeros(1,100) y(1,end-999:end)],'.-');
plot([y(2,1:1000) zeros(1,100) y(2,end-999:end)],'r.--');
legend('drift','no drift');
title('initial and last 1000 samples');

[P1, F1] = pwelch(y(1,:), hann(n), 0, n, Fs);
[P2, F2] = pwelch(y(2,:), hann(n), 0, n, Fs);
subplot(2,1,2); hold
plot(F1,10*log10(P1));
plot(F2,10*log10(P2),'r--');
legend('drift','no drift')
axis ([0 50 -20 100]);
grid
xlabel('MHz')

Phase continuity of NCO

Kadhiem Ayob December 9, 2011 Coded in Matlab
clear all; close all;

test = 1;                              %see comments below
n = 10000;                             %number of test samples

switch test
   case 1, z = linspace(.02,-.01,n);     %gentle zero crossing 
   case 2, z = linspace(.0005,-.0005,n); %excess dc at zero crossing
   case 3, z = randsrc(1,n,[.25 -.25]);  %random sweep, creates messy signal
   case 4, z = randsrc(1,n,[-.25:.001,.25]);  %random sweep, creates messy signal
   case 5, z = ones(1,n/4)*.01; z = [z -z z -z]; %discontinuity at zero crossing
end

%%%%%%%%%%%%%%%%%%%%%%%%% Finite NCO model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = 2^15 - 1;                          %max amplitude
M = 2^12;                              %NCO accumulator phase resolution
inc = round(M*z);                      %NCO accumulator phase increment
k = 2^12;                              %lut phase resolution (lut size), 
lsb = log2(M) - log2(k);               %LSBs discarded when addressing
lut = round(A*exp(j*2*pi*(0:k-1)/k));  %lut, one cycle cos/sin data

ptr = 0;
addr = 0;
for i = 1:n
    y(i) = lut(addr+1);                %add 1 for matlab LUT
    ptr = mod(ptr + inc(i), M);        %phase accumulator(ramp)
    addr = round(ptr/2^lsb);           %discard LSBs
    addr(addr >= k) = addr - k;        %check address overflow
end

%display output
plot(real(y),'-');hold
plot(imag(y),'r-');

I/Q and its conjugates

Kadhiem Ayob December 3, 2011 Coded in Matlab
clear all; close all;

fc = .019;  % frequency relative to Fs of 1(-.5 ~ +.5)
phase = 0;  % degrees

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
shift = floor(phase*1/(fc*360));
x1 = exp(j*2*pi*(0+shift:2^16-1+shift)*fc);   % I,Q
x2 = complex(real(x1),-imag(x1));             % I,-Q
x3 = complex(-real(x1),imag(x1));             %-I,Q
x4 = complex(imag(x1),real(x1));              % Q,I
x5 = complex(imag(x1),-real(x1));             % Q,-I
x6 = complex(-imag(x1),real(x1));             %-Q,I

f1 = fftshift(fft(x1,2^16));
f2 = fftshift(fft(x2,2^16));
f3 = fftshift(fft(x3,2^16));
f4 = fftshift(fft(x4,2^16));
f5 = fftshift(fft(x5,2^16));
f6 = fftshift(fft(x6,2^16));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = linspace(-.5,.5,2^16);
figure(1);
subplot(4,1,1);hold;
plot(f,20*log10(abs(f1)),'.-');
plot(f,20*log10(abs(f2)),'r-');
plot(f,20*log10(abs(f3)),'g--');
legend('I,Q','I,-Q','-I,Q')
ylabel('dB')

subplot(4,1,2);hold
plot(f,angle(f1),'.-');
plot(f,angle(f2),'r--');
plot(f,angle(f3),'g--');
legend('I,Q','I,-Q','-I,Q')
ylabel('rad')

subplot(4,1,3);hold;
plot(f,20*log10(abs(f4)));
plot(f,20*log10(abs(f5)),'r-');
plot(f,20*log10(abs(f6)),'g--');
legend('Q,I','Q,-I','-Q,I')
ylabel('dB')

subplot(4,1,4);hold
plot(f,angle(f4));
plot(f,angle(f5),'r--');
plot(f,angle(f6),'g--');
legend('Q,I','Q,-I','-Q,I')
ylabel('rad')

z = 1:ceil(1/abs(fc));
figure(2)
subplot(6,1,1);hold
plot(real(x1(z)),'.-');
plot(imag(x1(z)),'r.-');
legend('I','Q');

subplot(6,1,2);hold
plot(real(x2(z)),'.-');
plot(imag(x2(z)),'r.-');
legend('I','-Q');

subplot(6,1,3);hold
plot(real(x3(z)),'.-');
plot(imag(x3(z)),'r.-');
legend('-I','Q');

subplot(6,1,4);hold
plot(real(x4(z)),'.-');
plot(imag(x4(z)),'r.-');
legend('Q','I');

subplot(6,1,5);hold
plot(real(x5(z)),'.-');
plot(imag(x5(z)),'r.-');
legend('Q','-I');

subplot(6,1,6);hold
plot(real(x6(z)),'.-');
plot(imag(x6(z)),'r.-');
legend('-Q','I');