0

I have a signal with fc = 308 kHz and B = 30 kHz and some white noise added. The signal is sampled at samplFreq = 4fc.

chirp = {-0.0459056, -0.0473799, 0.0596458, 0.062499, -0.0756196, -0.0807912,
0.0937236, 0.102535, -0.113734, -0.127962, 0.135298, 0.15724, 
-0.157934, -0.190465, 0.181032, 0.227647, -0.203865, -0.268702, 
0.225606, 0.313447, -0.245347, -0.361593, 0.262129, 0.412745, 
-0.27497, -0.466406, 0.282897, 0.521979, -0.284988, -0.578774, 
0.2804, 0.63602, -0.26841, -0.692877, 0.24844, 0.748446, -0.220089, 
-0.801787, 0.183156, 0.851936, -0.137649, -0.897917, 0.0838031, 
0.938762, -0.0220736, -0.973525, -0.0468658, 1.0013, 0.122138, 
-1.02124, -0.202684, 1.03256, 0.287287, -1.03457, -0.374607, 1.02667, 
0.463214, -1.00838, -0.551628, 0.979336, 0.638359, -0.939322,
-0.721943, 0.888263, 0.800982, -0.826246, -0.874183, 0.75352, 
0.940383, -0.67051, -0.998577, 0.577813, 1.04794, -0.476203, 
-1.08786, 0.366628, 1.11789, -0.250201, -1.13782, 0.128191, 1.14765, 
-0.00200695, -1.14753, -0.126822, 1.13781, 0.256669, -1.119, 
-0.38584, 1.09173, 0.512603, -1.05673, -0.635223, 1.0148, 0.751999, 
-0.966813, -0.8613, 0.913627, 0.961605, -0.856113, -1.05154};

I use the code below to calculate the power spectrum

ts = If[EvenQ[Length[chirp]], timeSeries, Most[chirp]];
tslen = Length@ts;
power0 = PeriodogramArray[ts];

Now, how can I determine the corresponding frequency? I use this for positive frequency

freqPos = Table[(k samplFreq)/tslen, {k, Range[0, tslen/2. - 1]}];

If I just use freqNeg = - freqPos, I will have two zeros with different power values.

Mahdi Razaz
  • 337
  • 1
  • 5
  • Do you actually want the power spectrum or just the modulus of the spectrum? (This post)( https://mathematica.stackexchange.com/a/85167/12558) – Hugh Aug 26 '21 at 06:43

1 Answers1

1

Consider using Fourier instead of PeriodogramArray

You can use numpy definition for frequencies:

Clear[fftfreq] ;
fftfreq[n_Integer, d_Real:1.0] := 1/d * 1/n *  If[
    EvenQ[n],
    Join[Range[0, n/2-1], Range[-n/2, -1]],
    Join[Range[0, (n-1)/2-1], Range[-(n-1)/2, -1]]
] ;

Compare:

frequency = 0.1 ;
signal = Cos[2*Pi*frequency*Range[1000]] + 0.1*Cos[2*2*Pi*frequency*Range[1000]] + 0.05*Cos[3*2*Pi*frequency*Range[1000]] ;
range = fftfreq[1000] ;
fourier = Transpose[{range, Abs[Fourier[signal, FourierParameters -> {1, -1}]]}] ;
session = StartExternalSession["Python"] ;
function = ExternalFunction[
session,
"def function(signal):
  import numpy as np
  fourier = np.fft.fft(signal)
  range = np.fft.fftfreq(len(signal))
  return np.transpose([range, np.abs(fourier)])
"
] ;
npfourier = Normal[function[signal]] ;
ListPlot[{fourier, npfourier}, PlotRange->All, AspectRatio->1/2, PlotTheme->"Detailed", Filling->Axis, PlotStyle->{Directive[PointSize[Large],Black], Red}]
DeleteObject[session] ;

enter image description here

I.M.
  • 2,926
  • 1
  • 13
  • 18