Code Snippets Submitted by kaz
Even/Odd FIR filter structure
% odd/even filter structure verification
% two filter structures compared
% (1) direct use of filter function
% (2) even/odd split structure
x = randn(1,2048);
h = fir1(23,.3);
y1 = filter(h,1,x); % direct filter for reference
he = h(1:2:end); % even coeff, relative to hardware definition
ho = h(2:2:end); % odd coeffs
xe = x(1:2:end); % even samples
xo = x(2:2:end); % odd samples
f1 = filter(ho,1,xe); %subfilter 1
f2 = filter(he,1,xo); %subfilter 2
f3 = filter(he,1,xe); %subfilter 3
f4 = filter(ho,1,xo); %subfilter 4
%add up subfilters, advancing f3 one stage
ye = [0 f1] + [0 f2];
yo = [f3 0] + [0 f4];
%interleave even/odd into one stream
y2 = [y1 0 0];
y2(1:2:end) = ye;
y2(2:2:end) = yo;
%quantise and compare both ref and subfilters
y1 = round(y1*2^15);
y2 = round(y2*2^15);
plot(y1 - y2(2:end-1))Ideal interpolation filter
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');Phase drift of NCO
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
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-');Resampling filter performance
%%%%%%%%%%%%%%%%%% 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');Circular Convolution
%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')16QAM Modem model (Basic)
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')Reading text files 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])Power Computation of a Digital Stream
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);Vector alignment for verification
%"align_run" is top level for testing "align" function
function align_run
for i = 1:1
L1 = 3400; %length of first vector
L2 = 4023; %length of second vector
%generate two random vectors
x1 = randn(1,L1);
x2 = randn(1,L2);
%generate random index values for two identical segments
r1 = round(rand(1)*L1);
r2 = round(rand(1)*L2);
r1(r1 == 0) = 1;
r2(r2 == 0) = 1;
%length of identical segments as percent of vector length
n = round(50/100 * min([L1 L2]));
n(n >= min([L1 L2])) = min([L1 L2])-1;
while r1+n > L1, r1 = r1-1; end;
while r2+n > L2, r2 = r2-1; end;
%create identical segments at random location
x1(r1:r1+n) = x2(r2:r2+n);
%call align function
[y1,y2] = align(x1,x2);
%check alignment
subplot(2,1,1);
plot(y1);hold
plot(y2,'r--') ;
legend('y1','y2')
subplot(2,1,2);
plot(y1 - y2,'g')
legend('(y1 - y2)')
end
return
%%%%%%%%%%%%%%%%% align %%%%%%%%%%%%%%%%%%
function [y1,y2] = align(x1,x2)
%force a column orientation
x1 = x1(:);
x2 = x2(:);
L1 = length(x1);
L2 = length(x2);
L = max([L1 L2]);
%equalise vector lengths
if L1 > L2, x2 = [x2; zeros(L1-L2,1)]; end;
if L1 < L2, x1 = [x1; zeros(L2-L1,1)]; end;
%apply correlation and find index of maximum output
y = xcorr(x1,x2);
n = find(y == max(y)) ;
n = mod(n,L)+1;
%align x1 vector with x2
y1 = [x1(n:end); x1(1:n-1)] ;
y2 = x2;
returnCircular Convolution
%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
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
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
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
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
%%%%%%%%%%%%%%%%%% 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)
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
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
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
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');
