0 个投票

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')

请先登录，再回答此问题。

 采纳的回答

Paul
Paul 2022-1-24
编辑：Paul 2022-1-24

0 个投票

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).

6 个评论
显示 4更早的评论 隐藏 4更早的评论

Mohamed Ahmed
Mohamed Ahmed 2022-1-24
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
Paul 2022-1-25
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
Mohamed Ahmed 2022-1-25
编辑：Mohamed Ahmed 2022-1-25
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
Paul 2022-1-25
Ran in:
Ran in:
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
Mohamed Ahmed 2022-1-25
I'm trying to have a flat spectrum, which should be one of the properties of this signal
Paul
Paul 2022-1-25
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?

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by