0
$\begingroup$

I am translating my code from Python to Mathematica. Here is a piece of my Python code:

power_spectrum_k = np.zeros(n, float)
for k in range(1, int(n/2+1)):
    power_spectrum_k[k] = math.exp(-(2*math.pi*k*sigma/L)*(2*math.pi*k*sigma/L))

for k in range(int(n/2+1), n):
    power_spectrum_k[k] = power_spectrum_k[n-k]

inverse_transform2 = np.zeros(n, float)
inverse_transform2 = n*np.fft.ifft(power_spectrum_k)

CSI = zeros((n, n))
for i in range(n):
    for j in range(n):
        if abs(j-i) < (n/2):
            CSI[i][j] = inverse_transform2[abs(i-j)]
        else:
            CSI[i][j] = 0

How can I do that in Mathematica? In particular I don't know how to create the equivalent of the power_spectrum_k array, first, and then how to fill the CSI matrix.

Edit

This is what I have tried to do so far:

n = 64;
L = n;
sigma = 3;
nyquist = n/2 + 1;
sampling = 8;
mu = 0.0;
PowerSpectrum = Table[0, {n}];
PowerSpectrum[[0]] = 0.0;

For[i = 1, i <= (n/2), i++, 
  PowerSpectrum[i] = Exp[-(2*Pi*i*sigma/L)*(2*Pi*i*sigma/L)]];
For[i = (n/2 + 1) , i <= n, i++, 
  PowerSpectrum[i] = Exp[-(2*Pi*(n - i)*sigma/L)*(2*Pi*(n - i)*sigma/L)]];    
InverseTransformPowerSpectrum = InverseFourier[PowerSpectrum];
CSI = SparseArray[{i_, j_} /; 
  Abs[j - i] < (n/2) -> InverseTransformPowerSpectrum [[Abs[i - j]]], {n, n}]
$\endgroup$
6
  • $\begingroup$ You mean you want to translate power_spectrum_k = np.zeros(n, float) to Mathematica? I would use Table command, as in Table[0,{n}], for the whole thing, may be n=10; Table[Exp[-(2 Pi k w/L0)^2], {n}] ps. why not use ^2 instead of writing the same expression 2 times? $\endgroup$
    – Nasser
    Apr 12, 2015 at 5:53
  • $\begingroup$ All of your questions so far have been "Here's my Python code - how is this done in Mathematica..." with only one showing any effort to do this in Mathematica yourself. This is not a translation service. $\endgroup$
    – ciao
    Apr 12, 2015 at 5:53
  • 2
    $\begingroup$ You should not be translating code as is from Python to M. You should try to translate the algorithm itself if possible. Many operations in M can be done with one or 2 lines, instead of using For and so on. M is functional programming language, while Python is procedural build in an OO framework (like Java is). $\endgroup$
    – Nasser
    Apr 12, 2015 at 6:05
  • $\begingroup$ @Nasser I cannot do what you suggest in your first comment, because the values of the array are not always the same function of the index: they change from half of the array onwards. How can I deal with this issue? $\endgroup$
    – johnhenry
    Apr 12, 2015 at 6:09
  • $\begingroup$ Are we talking about this part: (2*math.pi*k*sigma/L)*(2*math.pi*k*sigma/L) ? Is this not the same as (2*math.pi*k*sigma/L)^2? This is what I mean. As for the table command index, yes, I see you have it going to int(n/2+1) and I used n so you just need to change the index. to one that corresponds to rounding n/2 and adding 1 to it. See help on Round, Ceiling, Floor etc.. $\endgroup$
    – Nasser
    Apr 12, 2015 at 6:15

1 Answer 1

2
$\begingroup$

As it's already mentioned in the comments, you need to define appropriate functions then use them in a functional form. For example, consider this function:

n = L= 64;
sigma = 3;
nyquist = n/2 + 1;
sampling = 8;
mu = 0.0; 
powerspectrum[i_] := 
  Piecewise[{
  {0, i == 0}, 
  {Exp[-(2*Pi*i*sigma/L)^2], 0 < i <= n/2}, 
  {Exp[-(2*Pi*(n - i)*sigma/L)^2], n/2 < i <= n}
  }]

If I have followed your equations correctly, the above function gives the value for any i from $0$ to $n$ for powerspectrum. Or even the above form can be written more neatly:

n = L= 64;
sigma = 3;
nyquist = n/2 + 1;
sampling = 8;
mu = 0.0;
g[exp_]:= Exp[-(2*Pi*exp*sigma/L)^2]
powerspectrum[i_] := Piecewise[{{0, i == 0}, {g[i], 0 < i <= n/2},{g[n - i], n/2 < i <= n}}]

Now you can use Table to make a list (array) of power values:

pts=Table[powerspectrum[i],{i,0,n}];
$\endgroup$
4
  • $\begingroup$ OK thanks this seems to work, but now, if I want to fill the CSI matrix according to the constraints that I have (see Python code) how can I do that in Mathematica? I guess a first step would be to perform the Inverse Fast Fourier Transform of what you called 'pts', which I tried to do using 'InverseFourier[pts]'...but then I don't know how to fill the matrix CSI following the constraints...Should I use Table again?If yes, how? $\endgroup$
    – johnhenry
    Apr 12, 2015 at 21:15
  • $\begingroup$ for example, I have tried with the following: $\endgroup$
    – johnhenry
    Apr 12, 2015 at 21:31
  • $\begingroup$ CSI = Table[0, {n}, {n}]; For[i = 1, i < n, ++i, For[j = 1, j < n, ++j, If[Abs[j - i] < (n/2), CSI[[i, j]] = inverse[[Abs[j - i]]], 0]]] $\endgroup$
    – johnhenry
    Apr 12, 2015 at 21:32
  • $\begingroup$ But I get the error: "Part specification inverse[[1]] is longer than depth of object." $\endgroup$
    – johnhenry
    Apr 12, 2015 at 21:32

Not the answer you're looking for? Browse other questions tagged or ask your own question.