How to plot Gauss sums ?

I'm trying to plot the Gauss sums according to the equation shown in the image s(t), but I keep receiving errors.

Can you please show me what am I doing wrong ?

clear all
close all
clc
%%
Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1024;                     % Length of signal
t = 2*(0:L-1)*T;               % Time vector
x = 0;

k = 0;
s = 0;
p = primes(L);

% s(t) = cumsum((k/p)(1:length(p)-1)).*exp(1i*k*t);

for k=1:p-1
    
s(t) = s(t) + (k/p).*exp(1i*k*t);
    
end

figure
subplot(2,2,1)
plot(t,s)
title('signal')

Answer 1

[Paul]

There are at least two problems with the code. Assuming you want to compute numerical values of s, the code can't reference s in a functional form, like s(t). For numerical objects, the "t" in s(t) is an index, which should be numerical integers or a logical. So for one value of p, you need to compute one value of an array s for each value of t, and then for each value of t, you need that inner loop over k. If you want more than one value of p, you'll need to make s a 2D array (one dimension for t and the other for p). Here is modification to the code for a single value of p

clear all
close all
clc
%%
Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1024;                     % Length of signal
t = 2*(0:L-1)*T;               % Time vector
x = 0;
k = 0;
s = 0;
p = primes(L);
p = 13; % odd prime
% s(t) = cumsum((k/p)(1:length(p)-1)).*exp(1i*k*t);
s = 0*t; % initialzize t
for jj = 1:numel(t)
for k=1:p-1
    
s(jj) = s(jj) + (k/p).*exp(1i*k*t(jj));
    
end
end

It's likely that at least the inner loop could be vectorized. However, this code has another problem. As stated in the text of the Question, the term (k/p) does not mean "k divided by p." Rather it is the Legendre symbol L(k,p). So you'll need a function that computes L(k,p).

Check out

doc jacobiSymbol

which would be perfect to use if p is an odd prime. But the problem statement isn't limited to only odd primes. However the wikipedia page for Legendre symbol seems to imply that it is only defined for odd primes. OTOH, the Mathworld page doesn't quite say that p is odd prime is assumed in the defintion (it might just be poorly worded).

[Mohamed Ahmed]

I tried using the code you provided and add the jacobiSymbol function, but I'm still stuck, as errors keep appearing to me.

Can you check if I did something wrong ?

clear all
close all
clc
%%
Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1024;                     % Length of signal
t = 2*(0:L-1)*T;              % Time vector
x = 0;
k = 0;
s = 0;
p = primes(L);
s = 0*t; % initialize t

for jj = 1:numel(t)
for k=1:p-1

J = jacobiSymbol(k,p);
m = powermod(k,(p-1)/2,p);
checkPrime = all(mod(J,p) == m)
isprime(p)
s(jj) = s(jj) + J.*exp(1i*k*t(jj));
    
end
end

figure
subplot(2,2,1)
plot(t,s)
title('signal')

[Paul]

The code is set up to only deal with one value of p at a time. But your code is computing 172 values of p, and then the "for k = 1:p-1" results in something weird, and some of those values of p are even, but jacobiSymbol only allows for p to be odd. Here is the same code I posted above, but for one value of p and using jacobiSymbol

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1024;                     % Length of signal
t = 2*(0:L-1)*T;              % Time vector
x = 0;
k = 0;
s = 0;
%p = primes(L);
s = 0*t; % preallocate s
p = 13;  % single value of p!
for jj = 1:numel(t)
for k=1:p-1
J = jacobiSymbol(k,p);
% m = powermod(k,(p-1)/2,p);
% checkPrime = all(mod(J,p) == m)
% isprime(p)
s(jj) = s(jj) + J.*exp(1i*k*t(jj));
    
end
end

This code can probably be much improved, but first it's important for you to get it working.

If you want to compute s for multiple, odd prime values of p, you can add a second dimension to s and another loop over the values of p

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1024;                     % Length of signal
t = 2*(0:L-1)*T;              % Time vector
x = 0;
k = 0;
%s = 0;
%p = primes(L);
pvalues = [11 13];  % odd and prime!
s = zeros(numel(t),numel(pvalues));
for np = numel(pvalues)
    p = pvalues(np);
    for jj = 1:numel(t)
        for k=1:p-1
            J = jacobiSymbol(k,p);
            s(jj,np) = s(jj,np) + J.*exp(1i*k*t(jj));
        end
    end
end

[Mohamed Ahmed]

I added more prime numbers and removed the for loop "for jj = 1:numel(t)"

but now the plot is not what is expected from the magnitude response which should be flat.

clear all
close all
clc
%%
Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1024;                     % Length of signal
t = 2*(0:L-1)*T;              % Time vector
x = 0;
k = 0;
%s = 0;
%p = primes(L);
pvalues = [2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73];  % odd and prime!
% 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 127
s = 0;
for np = numel(pvalues)
    p = pvalues(np);
        for k=1:p-1
            J = jacobiSymbol(k,p);
            s = s + J.*exp(1i*k*t);
        end
end

figure
subplot(2,2,1)
plot(t,s)
title('signal')


%%
X = x;
n = 2^nextpow2(L);
f=0:(Fs/n):(Fs/2-Fs/n);

VxModel = s;
tRangeModel = t;
VxModel = VxModel - mean(VxModel);           %When calculating the FFT and divide by number of samples, indeed you get the mean value of your signal. In order to remove that, please apply the following before the FFT calculation.

VFFTModel = fft(VxModel,n);
VFFTModel = transpose(VFFTModel(:));
NVFFT   = length(VFFTModel);
MGVFFT2 = abs(VFFTModel/NVFFT);
MGVFFT1 = 2*MGVFFT2(1:floor(NVFFT/2));

PHVFFT2 = angle(VFFTModel);
PHVFFT1 = PHVFFT2(1:floor(NVFFT/2));

MGVFFT1Model = MGVFFT1;%(FListIndexModel);
PHVFFT1Model = PHVFFT1;%(FListIndexModel);
VxFFTModel = (MGVFFT1Model.*exp(j*PHVFFT1Model));

subplot(2,2,2)
semilogx(f,MGVFFT1Model)
title('Magnitude')
subplot(2,2,3)
semilogx(f,rad2deg(PHVFFT1Model))
title('Phase')

VRecalFFT = [(MGVFFT1Model*NVFFT/2.*exp(j*PHVFFT1Model))  0 ...
    (MGVFFT1Model(end:-1:2)*NVFFT/2.*exp(-1*j*PHVFFT1Model(end:-1:2)))];% zeros(1,NVFFT/2)]; 
VxReModel1 = ifft(VRecalFFT,'symmetric');
VxReModel2 = ifft(VFFTModel);

subplot(2,2,4)
plot(tRangeModel,VxModel,'b.',...
    tRangeModel,VxReModel1,'r.',...
    tRangeModel,VxReModel2,'b','linestyle','-.')

title('Reconstructed signal')

[Paul]

My apologies, the code I posted above had an error. Here is the corrected code:

Fs = 1000;                    % Sampling frequency
T = 1/Fs;                     % Sampling period
L = 1024;                     % Length of signal
t = 2*(0:L-1)*T;              % Time vector
x = 0;
k = 0;
pvalues = [11 13];  % odd and prime!
s = zeros(numel(t),numel(pvalues));
for np = 1:numel(pvalues) % corrected line
    p = pvalues(np);
    for jj = 1:numel(t)
        for k=1:p-1
            J = jacobiSymbol(k,p);
            s(jj,np) = s(jj,np) + J.*exp(1i*k*t(jj));
        end
    end
end

Now s is a two dimensional matrix where the the rows correspond to the values of t and the columns to each value of t. s is complex. Is that what you expected?

whos pvalues t s

  Name            Size              Bytes  Class     Attributes

  pvalues         1x2                  16  double              
  s            1024x2               32768  double    complex   
  t               1x1024             8192  double              

Make a plot

figure
plot(t,abs(s))

Now at this point I think you should be able to operate on each column of s in whatever fashion you need. If you can't get the the next steps to work, please explain what you're trying to do with s before changing any of the code in this comment.

Also, p = 2 is not odd.

[Mohamed Ahmed]

I'm trying to have a flat spectrum, which should be one of the properties of this signal

[Paul]

s(t) is a sum of complex exponentials. The Fourier transform of a complex exponential is Dirac delta, so the spectrum of s(t) should be an train of p-1 delta functions, shouldn't it?