Phase Error Growth in Digital Mixer Signal

Question

I am working on a code where digital mixing is to be performed for down-conversion of incoming signal. The mixer signal is created using the phase accumulation approach. However it was observed that the phase error between a signal created with phase accumulation vs that created with a fixed duration time signal had some error which was observed to grow. Can anyone explain the difference and suggest how to avoid this. Code and example result below.

clear all;close all;
fs=25e3;
sim_time=40e-3;
f_mix=5e3;
t=0:1/fs:100*sim_time-(1/fs);
sig_mix_lref=cos(2*pi*f_mix*t)+j*sin(2*pi*f_mix*t);%create mixer reference
ph_ref=wrapToPi(2*pi*f_mix*t);%phase ref for comparison with accumulated phase

ph_acc(1)=0;%init phase
ph_inc=2*pi*f_mix*1/fs;%phase increment at desired rate
ind=2;

for i=1:2048:length(sig_mix_lref)-2048
    for k=1:2048
        sig_mix_l_c(k)=exp(1*j*ph_acc(ind-1));%create mixer block later appended for comparison with reference
        ph_acc(ind)=ph_acc(ind-1)+ph_inc;
        ph_acc(ind)=wrapToPi(ph_acc(ind));
        %do block mixing here with incoming data stream
        ind=ind+1;      
    end
    sig_mix_l(i:i+2048-1)=sig_mix_l_c;
end

sig_mix_lref=sig_mix_lref(1:length(sig_mix_l));

ph_ref=ph_ref(1:length(ph_acc));

figure;plot(real(sig_mix_lref)-real(sig_mix_l),'-*r')
title('mixer comparison error real part')

figure;plot(imag(sig_mix_lref)-imag(sig_mix_l),'-*b')
title('mixer comparison error imag part')

figure;plot(ph_ref-ph_acc,'-*r')
title('mixer accum phase error wrt to ref phase') 

Although the error scale is negligible, but can anyone explain the reason for growth, how to resolve it?

Answer 1

This is numerical noise. You are using floating point math which is never "exact" but only approximate. The problematic line here is

ph_ref=wrapToPi(2*pi*f_mix*t);%phase ref for comparison with accumulated phase

The argument into your phase wrapping function can get large and the larger the argument, the larger the phase wrapping error will be.

You can simulate the effect by taking any angle add a large number of multiples of $2\pi$ unwrap it and look at the difference to the original angle. The result looks roughly like the graph below (for 64 bit double precision floating point).

Note that the phase accumulation is better. It always keeps the phase between $0$ and $2\pi$ and the unwrapping can simply be done with if phase > 2*pi, phase = phase - 2*pi; end which avoids division and/or modulo calculations.

Here is the code to generate this graph

%% look at phase wrapping error as a function of length

a = 2*pi*(0.5:999.5)'/1000; % test vector with 100 angles
numberOfPeriods = 10.^(0:15)'; % number of peroids to unwrap
n = length(numberOfPeriods);
errorDb = zeros(n,1);
% loop over all prders of magnitudes
for i = 1:n
  x = numberOfPeriods(i)*2*pi;
  % wrap angle to [0, 2*pi|
  wrappedAngle = (x+a)-2*pi*floor((x+a)./(2*pi));
  % calculate  angle difference
  angleDiff = wrappedAngle-a;
  % relative error in dB
  errorDb(i) = 10*log10(mean(angleDiff.^2)./mean(a.^2));
  
end

%% and plot the result
clf;
h = semilogx(numberOfPeriods,errorDb);
set(h,'Linewidth',2);
grid('on');
xlabel('Number of Periods in vector');
ylabel('Relative error in dB');
title('Phase wrapping error');