Code covered by the BSD License  

Highlights from




15 Jan 2014 (Updated )

Phase-locked loop.m this is an mfile to calculate the same outputs of the PLL simulink block.

function [f,f_time,f_interp,time_Zeros,Wt,Sine,Cose,Signal_Zeros,Signal_Limits,PLL_time] = PLL_final_cal(Signal,time,tol,option)
% This function will locate the zero crossing in your input signal and
% calculate the frequency of that input signal.
% This function will also generate your input signal 'Wt' function and
% regenerate an pu sin & cos signal just like an PLL_Block in simulink.
% the needed inputs:
% 1- Signal: is your input y-axi data (in case of sine wave those are the
% Amplitude)
% 2- time: is your input x-axi data which is the time space of the input
% signal
% tol: is the tolerance in which zeros points will be selected. and it
% should be given in Hz.
% for example if tol = 1 [Hz] is means any Zeros(i+1) point which distance
% from the point before Zeros(i) is in this tolerance band of
% [1/(fn-tol),1/fn,1/(fn+tol)] in S will be accepted as the next frequency
% of f(i+1).
% the available outputs:
% 1- f: is the frequency of the input Signal calculated at every Zeros(i)
% points which means f length is the same as the number of Zeros found in
% the input Signal divided by 2 (in the Sine wave the half wave part is
% neglected). mins 1 (due to the "diff" function effect).
% 2- f_time: is the time signals which can be used to "plot" this f signal
% and returns the times where the Zeros has occured and the frequency were
% calculated.
% f_interp & time_Zeros: are the interpulation of the previous signals in
% order to get the output frequency with the same size of the input Signal.
% 3- Wt, Sine & Cose: are the same output of the PLL function block in
% simulink. which means the are the pu represntation of the input Signal
% where there is no harmonics. yet have the same frequency of the input
% Signal.
% 4- Signal_Zeros: is the cut version of the input Signal (to synchronize
% on the first or second Zero crossing)
% 5- Signal_Limits: is the indx representation of the input signals with
% only two values [Start_Value End_Value], in which the complete output was
% created.
% 6- option: use this variable to do the calculation accordirng to the PLL
% simulink block (When option = 1) or to build up your own Wt data from the
% zero crossing points (When option = 2)
% this function file will need the following function files in order to
% work:
% 1- Fast_Harmonics.m: is one is used to extract the fundamental frequency
% out of the input signal.
% 2-  Zeros_finding_Rowdata_final.m: this function will be used to extract
% the zero crossing of the input signal and then internally takes only the
% zeros points which fulfill the 'tol' argument mentioned before. so for
% example if you have an input signal where too much harmonics are inside
% and some missleading zero crossing in there is function will try to
% remove those unwanted zeros points.
% Author: Aubai Al Khatib datum: 15.01.2014
Delta_Time_Orginal = time(2)- time(1);
F_Sample = abs(1/Delta_Time_Orginal);
%---------------------- The Harmonics Spectrum of the complet input signal ----------------------------
Tp = length(time);
T_S = abs(Tp*Delta_Time_Orginal);
[nyquist_number_I1,Bezugsfrequenz_I1,Signal_amp_I1] = Fast_Harmonics (Signal,F_Sample,T_S);
%plot([0:nyquist_number_I1]*Bezugsfrequenz_I1,Signal_amp_I1);ylabel('amplitude Currents [A]');title('Harmonic spectrom');hold on
F_Spec = [0:nyquist_number_I1]*Bezugsfrequenz_I1; %#ok<*NBRAK>
Fn = round(F_Spec(Signal_amp_I1 == max(Signal_amp_I1)));
[~,Zeros_f,Zeros_f_indx_final] = Zeros_finding_Rowdata_final(Signal,time,Fn,tol);
%[Zeros,k] = Zeros_finding(Signal,time);
Zeros = Zeros_f;
k = Zeros_f_indx_final;
f = 1./(2*diff(Zeros));
f_time = Zeros(1:length(f));
Test_indx = k(1)+2;
if Zeros(1) < 1/(Fn*2) && Signal(Test_indx) < 0
    Signal_Zeros = Signal(k(2):k(end-1));
    time_Zeros = time(k(2):k(end-1));
    Signal_Zeros = Signal(k(1):k(end-1));
    time_Zeros = time(k(1):k(end-1));
f_interp = cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), f, diff(k),'uniformoutput',0));
time_f_interp = time(k(1):k(end)-1);
time_Zeros_1 = time_Zeros;
time_Zeros = time_f_interp;
if option
    W = 2*pi*f_interp;
    Wt = cumtrapzt(time_Zeros,W);
    Wt = mod(Wt,2*pi);
    if Zeros(1) < 1/(Fn*2) && Signal(Test_indx) < 0
        Start_Indx = k(2);
        End_Indx = k(end-1);
        Sine = -sin(Wt);
    Cose = -cos(Wt);
        Start_Indx = k(1);
        End_Indx = k(end-1);
        Sine = sin(Wt);
    Cose = cos(Wt);
    Signal_Limits = [Start_Indx, End_Indx];
    PLL_time = time_f_interp;
    Zeros_one_period = Zeros(1:2:end);
    D_Zeros = kron(Zeros_one_period,[1 1])';
    Test = 2:2:length(D_Zeros);
    D_Zeros(Test) = D_Zeros(Test)+0.00000001;
    A_max = 0;
    A_min = 2*pi;
    A = [A_min; A_max];%*(pi/180);
    R_Z_G = repmat(A,length(D_Zeros),1);
    R_Z_G_1 = R_Z_G(1:length(D_Zeros));
    %R_Z_G_2 = R_Z_G(1:length(Zeros));
    R_Z_G_interp = interp1(D_Zeros,R_Z_G_1,time,'linear');
    %R_Z_G_interp = interp1(R_Z_G_1,time,'linear');
    R_Z_G_interp_test = R_Z_G_interp(~isnan(R_Z_G_interp));
    Wt = R_Z_G_interp_test;
    Test_indx = k(1)+2;
    if Zeros(1) < 1/(Fn*2) && Signal(Test_indx) < 0
        Sine = -sin(R_Z_G_interp_test);
        Cose = -cos(R_Z_G_interp_test);
        Sine = sin(R_Z_G_interp_test);
        Cose = cos(R_Z_G_interp_test);
    Time_test = time(isnan(R_Z_G_interp) == 0);
    PLL_time = Time_test;
    Start_Indx = find(time == Time_test(1));
    End_Indx = find(time == Time_test(end));
    Signal_Limits = [Start_Indx, End_Indx];
    Signal_Zeros = Signal(isnan(R_Z_G_interp) == 0);
plots = 1;
if plots
    figure(); %#ok<*UNRCH>
    x_min = 200;
    x_max = x_min + 0.1;%time(end);%0.1;
    [AX,ha2,ha3] = plotyy(time,Signal,time_Zeros,Sine);grid on;  %#ok<*NASGU,*ASGLU>
    set(AX(1),'xlim',[x_min x_max])
    set(AX(2),'xlim',[x_min x_max])
    set(AX(1),'ylim',[-max(Signal) max(Signal)])
    plot(time_Zeros,f_interp);grid on;xlim([x_min x_max]);
    plot(time_Zeros,Wt);grid on;xlim([x_min x_max]);
    plot(time_Zeros,Sine);grid on;xlim([x_min x_max]);hold on
    plot(time_Zeros,Cose,'r');hold off

Contact us