How to plot periodogram from DSP.jl?

General Usage
, , ,
1

I am trying to plot a periodogram for the following data. However, I am not sure how can i plot a periodogram object ?

Code:

#= DataFrame
10×10 DataFrame
 Row │ Date        Open     High     Low      Close    Volume     AdjustedOpen  AdjustedHigh  AdjustedLow  AdjustedClose 
     │ Date        Float64  Float64  Float64  Float64  Float64    Float64       Float64       Float64      Float64       
─────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ 2020-08-21    41.3   41.34     40.71     41.01  2.55177e7       38.5015       38.5388      37.9515        38.2312
   2 │ 2020-08-24    41.28  42.5078   41.09     42.22  1.82089e7       38.4829       39.6275      38.3058        39.3592
   3 │ 2020-08-25    41.66  42.065    40.7      40.88  3.12856e7       38.8371       39.2147      37.9422        38.11
   4 │ 2020-08-26    40.64  40.835    40.0      40.01  2.6795e7        37.8862       38.068       37.2896        37.2989
   5 │ 2020-08-27    40.0   40.3      39.305    39.74  3.20343e7       37.2896       37.5693      36.6417        37.0472
   6 │ 2020-08-28    39.84  40.92     39.76     40.69  3.34771e7       37.1405       38.1473      37.0659        37.9329
   7 │ 2020-08-31    40.64  40.72     39.895    39.94  2.52596e7       37.8863       37.9608      37.1917        37.2337
   8 │ 2020-09-01    39.75  39.75     39.045    39.43  2.24782e7       37.0565       37.0565      36.3993        36.7582
   9 │ 2020-09-02    39.23  39.715    38.95     39.19  2.64134e7       36.5718       37.0239      36.3108        36.5345
  10 │ 2020-09-03    39.2   40.03     38.88     39.11  2.88172e7       36.5438       37.3176      36.2455        36.4599
=#
using DSP, DataFrames, Plots
data = df[:, :AdjustedClose]
a = periodogram(data)

Thanks look forward to the suggestions.

2

LPVSpectral.jl defines plot recipes for periodograms and spectrograms

julia> using DSP, LPVSpectral

julia> plot(periodogram(randn(1000)))

per
per1100×800 122 KB

4 Likes
3

@mdsa3d, one difficulty with your input data is that the time axis has a non-uniform sampling rate, which is a requirement of the FFTW. I.e., the Date column in the dataframe is not sampled at a regular time interval.

Strictly speaking, NFFT.jl or equivalent, should be used but AFAIK, there are no practical examples posted.

Code below using NFFT.jl was edited as well as example in post further down comparing NFFT.jl with FINUFFT.jl:

using Dates, DataFrames, NFFT, Statistics, Plots; gr()

colname = [:Date, :Open, :High, :Low, :Close, :Volume, :AdjustedOpen, :AdjustedHigh, :AdjustedLow, :AdjustedClose]

M =["2020-08-21"  41.3   41.34    40.71   41.01  2.55177e7  38.5015  38.5388  37.9515  38.2312;
     "2020-08-24"  41.28  42.5078  41.09   42.22  1.82089e7  38.4829  39.6275  38.3058  39.3592;
     "2020-08-25"  41.66  42.065   40.7    40.88  3.12856e7  38.8371  39.2147  37.9422  38.1100;
     "2020-08-26"  40.64  40.835   40.0    40.01  2.6795e7   37.8862  38.068   37.2896  37.2989;
     "2020-08-27"  40.0   40.3     39.305  39.74  3.20343e7  37.2896  37.5693  36.6417  37.0472;
     "2020-08-28"  39.84  40.92    39.76   40.69  3.34771e7  37.1405  38.1473  37.0659  37.9329;
     "2020-08-31"  40.64  40.72    39.895  39.94  2.52596e7  37.8863  37.9608  37.1917  37.2337;
     "2020-09-01"  39.75  39.75    39.045  39.43  2.24782e7  37.0565  37.0565  36.3993  36.7582;
     "2020-09-02"  39.23  39.715   38.95   39.19  2.64134e7  36.5718  37.0239  36.3108  36.5345;
     "2020-09-03"  39.2   40.03    38.88   39.11  2.88172e7  36.5438  37.3176  36.2455  36.4599]

df = DataFrame(M, colname)
df.Date = Date.(df.Date)
time = (df.Date .- df.Date[1]) ./ Day(1)   # non-uniform days vector
data = Complex.(df.AdjustedClose)
md = mean(data)
n = size(df,1)

# normalize time span of 14 days over [-0.5, 0.5)
tnrm = range(-0.5,0.5,n+1)[1:n]
p = plan_nfft(tnrm, n)
g = nfft(p, data .- md)
freq = (-n÷2:n÷2-1)/time[end]

p1 = plot(time, Real.(data), xlabel="Day", ylabel="AdjustedClose")
p2 = plot(freq[n÷2+1:n], abs.(g)[n÷2+1:n], xlabel="1/Day", ylabel="Spectral amplitude",
           title="using NFFT.jl")
plot(p1,p2,legend=false, markershape=:circle, ms=2)

NFFT_Type1

NB: depending on the actual problem, time interpolation + FFTW.jl might be enough.

1 Like
4

FWIW, please find below a practical example using FINUFFT.jl, another Julia package for NFFT computations, with results compared to NFFT.jl.

The code uses Dates and days to be in line with OP’s problem of an irregularly sampled time series.

using Dates, FINUFFT, Statistics, Random, Plots; gr()

Random.seed!(111)

# INPUT DATA:
Tp = 365                                          # period domain
t = sort(Day.(rand(1:2Tp, 200)))  ./ Day(1)       # non-uniform days vector
st = 40 .+ sin.(2π*(3/Tp)*t) + cos.(2π*(7/Tp)*t)  # data values
n = length(t)                                     # number of points

# DATA PRECONDITION:
tn = range(-π,π,n+1)[1:n]     # normalize to [-π,π)
st = Complex.(st)             # complex data vector type required
M  = 4*n                      # number of targets to compute the Fourier Transform

# COMPUTE (fast nufft of type-1):
tol = 1e-10
Ft = nufft1d1(collect(tn), st .- mean(st), +1, tol, M)
freq = (-M÷2:M÷2-1)/t[end]

# PLOT:
# The frequency peaks at 3/Tp (=0.008) and at 7/Tp (=0.019) are retrieved correctly:
p1 = plot(t, real(st), xlabel="Time (days)", ylabel="signal amplitude");
p2 = plot(freq[M÷2:M], abs.(Ft)[M÷2:M], xlabel="Frequency (1/day)", ylabel="Spectral amplitude",
          xlims=(0,0.05), title="using FINUFFT.jl");
plot(p1,p2,legend=false, markershape=:circle, ms=2)


# Compare to NFFT.jl. normalize time over [-0.5, 0.5):
# COMPUTE (fast nfft of type-1):
using NFFT
tn = range(-0.5,0.5,n+1)[1:n]
p = plan_nfft(tn, n)
Ft2 = nfft(p, st .- mean(st))
freq = (-n÷2:n÷2-1)/t[end]

# PLOT:
p1 = plot(t, Real.(st), xlabel="Day", ylabel="AdjustedClose")
p2 = plot(freq[n÷2+1:n], abs.(Ft2)[n÷2+1:n], xlabel="Frequency (1/day)", ylabel="Spectral amplitude",
          xlims=(0,0.05), title="using NFFT.jl")
plot(p1,p2,legend=false, markershape=:circle, ms=2)

The frequency peaks at 3/Tp (=0.008) and at 7/Tp (=0.019) are retrieved:

FINUFFT_Type1_irregular_time_sampling

NFFT_Type1_irregular_time_sampling

3 Likes
5

Thank you for this amazing suggestion, highly appreciate it!!!
It has cleared numerous doubts for me !!!

6

Fyi, the code above was updated and in one of the posts a comparison is provided between NFFT.jl and FINUFFT.jl results.

1 Like