code updated, now works

main.m

@@ -23,7 +23,7 @@
 %% Generation of synthetic signal
 A           = 1;        % mW/s       Amplitude of the synthetic signal
 cno         = 60;       % dbHz       C/N0 of the synthetic signal
-phi         = deg2rad(10);
+phi         = deg2rad(30);
 % phi = linspace(0, 2*pi, length(tVec));
 % phi = sin(2*pi*5*tVec);
 signal      = syntheticSignal(prn, kDelay, A, cno, tC, fIF, phi, tVec, fDoppler, B);
@@ -42,21 +42,23 @@
 
 %% PLL
 nSamples = tI*fs;
-theta = zeros(1,length(tVec));
+theta_hat = zeros(duration/tR,1);
 Vd = zeros(duration/tR,1);
 Vc = zeros(duration/tR,1);
 I = zeros(duration/tR,1);
 Q = zeros(duration/tR,1);
 k = 2;
 % iTs = 0:1/fs:tI-1/fs;
-iTs = 1/fs:1/fs:tI;
+iTs = (1/fs:1/fs:tI) - tI/2;
 for t=1:nSamples:length(signal)-nSamples
     In = signal(t:t+nSamples-1)     ...
         .*cm(t:t+nSamples-1)        ... %add doppler in code
-        .*cos(2*pi*(fIF+fDoppler)*tVec(t:t+nSamples-1) - theta(t:t+nSamples-1));
+        .*cos(2*pi*(fIF+fDoppler)*tVec(t:t+nSamples-1) ...
+        - 2*pi*Vc(k-1)*iTs - theta_hat(k-1));
     Qn = -signal(t:t+nSamples-1)    ...
         .*cm(t:t+nSamples-1)        ...
-        .*sin(2*pi*(fIF+fDoppler)*tVec(t:t+nSamples-1) - theta(t:t+nSamples-1));
+        .*sin(2*pi*(fIF+fDoppler)*tVec(t:t+nSamples-1) ...
+        - 2*pi*Vc(k-1)*iTs - theta_hat(k-1));
 
     I(k) = (1/nSamples)*sum(In);
     Q(k) = (1/nSamples)*sum(Qn); 
@@ -65,7 +67,7 @@
 
     switch order
         case 1 % First order
-            theta(t+nSamples:t+2*nSamples-1) = theta(t:t+nSamples-1) + Kpll * Vd(k);
+            theta_hat(k) = theta_hat(k-1) + Kpll * Vd(k);
         case 2 % Second order
             K1 = 60/23 * Bl * tI; 
             K2 = 4/9 * K1^2;
@@ -79,7 +81,8 @@
                     (K1+K2+K3) * Vd(k)/(2*pi*tI)  - ...
                     (2*K1+K2) * Vd(k-1)/(2*pi*tI) + ...
                     K1 * Vd2/(2*pi*tI);
-
+
+            theta_hat(k) = theta_hat(k-1) + 2*pi*Vc(k)*tI;
 %             theta(t+nSamples:t+2*nSamples-1) = theta(t:t+nSamples-1) + 2*pi*Vc(k)*iTs;
 %             avgTheta(k) = mean(theta(t+nSamples:t+2*nSamples-1))
         otherwise
@@ -106,7 +109,7 @@
 ylabel('Amplitude');
 title('Quantized signal');
 
-figure; plot(theta(nSamples/2:nSamples:end)*180/pi);
+figure; plot(theta_hat*180/pi);
 xlabel('Time (ms)'); ylabel('Phase (deg)');
 title('Carrier Phase Tracking');
 
@@ -135,10 +138,13 @@
 xlabel('Vd (deg)'); ylabel('Counts');
 title('Distribution of Vd');
 
-% TODO: Add qqplot
+figure;
+qqplot(rad2deg(VdConv))
+xlabel('Standard Normal Quantiles'); ylabel('Quantiles of Vd');
+title('QQ Plot of Vd vs Standard Normal');
 
 [Swelch, fwelch] = pwelch(VdConv, 512, 0, 512, 1/tI);
-figure
+figure;
 plot(fwelch./1e3, 10*log10(Swelch));
 xlabel('f (kHz)'); ylabel('PSD (dB/Hz)');
 title('Welch periodogram of the discriminator output');