Suggested algorithm for computing the magnitude spectrum and time-domain signal for the Complx FM equation

Hi, I am having some trouble getting some code to work properly and return the waveform of a complex FM signal and also return the spectrum computed using the FFT. My code thus far is as follows (notice: not all of the code works correctly) I am really stuck with this am would appreciate any prompt solution/suggestions.

clear all; %initialise parameters as used in paper I = [1 0.7 0.2]; omegac = 2*pi*100; theta = 0; phi = [0 0 0]; omegam = omegac;

K = length(I); % Number of Indexs T = 1; % Initialise period (T) TotalNumFreq = 0; % Initialise TotalUniqueFreq = 0; % Initialise

for index0 = 1:K % For the number of modulation indexes ki = ceil(1.2185*I(index0) +5.625) % equation 10 (see paper) Orders for bessel functions T = T*2*(ki-1) % Equations 12-15 (see paper) Waveform Period

    kivec(index0) = 2*ki - 1 % Evolution of index values
    TotalUniqueFreq = 1 + 2*sum(index0)*ki % Equation 20 (see paper) Total Number of Unique Frequencies

    end % End 1st for statement

% kf is a concatenation with kivec kf = [1 kivec] %kf is indexed from 1 to k+1

xi = zeros(T,K); for index1 = 1:T % For period (1 to T) for index2 = 1:K % For number of indexes tau(index2) = (index2-1)/T;

        if index1 < K-1 % If the value of T is less than K-1

            ind = cumprod(kf(1:index2));
xi(index1, index2)= 0.5(kivec(index2))*sawtooth(((2*pi)/ind(index2))*tau(index2)) - 0.5*sawtooth(((2*pi)/ind(index2))*tau(index2));

        else
            ind2 = cumprod(kf(1:K));
            xi(K, index2) = 0.5*(kivec(K))*sawtooth(((2*pi)/ind2(K))*tau(K)) + 0.5;

        end % End if statement

J(index1,index2) = besselj(xi(index1,index2),I(index2));

omega (index1, index2) = omegam * index2 * xi(index1, index2);

phi1 (index1, index2) = phi(index2) * xi(index1, index2);

    end % temporary end for testing

% end % End for statement

omegacap = cumsum(omega,2); % Frequencies Array

gamma = cumsum(phi1,2); % Phase Array

A = cumprod(J(1:index2),2);

end % temporary end for testing

TotalNumFreq = T; % Total number of frequencies is now T

for index3 = 1:TotalUniqueFreq % For the total number of unique frequencies

    FreqIndex = index3 - TotalUniqueFreq/2; % Find indexes for frequency
    indexF = []; % Add element to frequencies array

    for index4 = 1:TotalNumFreq % For the total number of frequencies

if omegacap(index4) == FreqIndex*omegam

           [indexF] = [indexF index4]; 

       end % End if
    end % End for statement

D(omegac + FreqIndex*omegam) = cumsum(A, indexF);

end % End for statement

 if resynthesis == True % If resynthesising

     d(t) = symsum(A*sin((omegac + omegacap)*t + gamma + theta), index1);

end % End if statement

Thanks in advance for any help here. I am happy to provide any additional information that might be required.

Richard.