I would like to know, how may I filter or separate this complex function into its simpler sub-functions (something like fourier filter).
I have a noisy dataset:
x = 1:1:10
y = [-1,5,0,2,-7,1,6,9,-2,3]
If I would like to find all simple functions which are combined to form this complex data.
What would be the correct way to execute this operation. Please suggest an approach to this probelm.
Thanks in advance!
For the Fourier decomposition of a signal as a sum of sinusoidal functions, you may take a look at example here.
To process your signal in the frequency domain, you can use the DSP.jl package to filter out the higher frequency components of your “noisy” dataset:
using DSP, Plots
x = 1:1:10
y = [-1.0,5,0,2,-7,1,6,9,-2,3.0]
fs = 1.0; # sampling frequency
fc = 0.35; # frequency cutoff, less than Nyquist
responsetype = Lowpass(fc, fs=fs)
designmethod = FIRWindow(hamming(10))
filty = filtfilt(digitalfilter(responsetype, designmethod), y)
plot(x, y, label = "input", legend=:topleft)
plot!(x, filty, label = "filtered")
PS: it goes without saying, the input array provided is pretty short.
Thanks @rafael.guerra for the explanation and link suggestions.
My goal was to convert a timeseries complex pattern into many simpler pattern.
As I am working with this dataset which has lot of variations. Hence it is computationally expensive to work with this complex pattern. Thus if i can somehow convert this complex pattern into number of simpler, smoother and small patterns then it wont be heavy on my system.
I have tried your fourier filter example and will start working from it but if you have some suggestion on how this can be achieved please do let me know.
Thanks for you help!
Can you tell us more about what type of data your are using ? Depending on that, you can choose a basis of analysis (or bank of filters) made of simple functions that you can use to filter your signal. The basis you choose will depend on the analytic properties of your signal (like is it varying fast or not) and what type of information you want to highlight.
Let’s note this basis \{\psi_i\}_i, 1\leq i \leq L. If we write y your data signal getting value y(t) at time t, by assuming that you want to keep the analysis in the time domain, you can have a new signal x_i(t) in the time domain by convoluting it with the i-th element of your basis:
Where * is the convolutional operator. The computation is simply the inverse Fourier Transform of the product of the Fourier Transforms \mathcal{FT}(y) and \mathcal{FT}(\psi_i) of y and \psi_i. This computation is fast this way using the Fast Fourier Transform (you want to use FFTW.jl in Julia).
Depending on what you chose for \psi_i you can compress x_i such that it has smaller size. I assume that your signal is real valued (\in \mathbb{R}) and not complex valued (\notin \mathbb{C}). Without getting into too much details, you can get a compression factor of n/(n-2*n_{\mathrm{zeros}}) where n is the size of your signal and n_{\mathrm{zeros}} the number of values close to zero the Fourier Transform \mathcal{FT}(\psi_i) of your function \psi_i have (outside a compact Fourier domain). Just to show you that to answer your problem we need some context on the type of data (is it coming from physical processes, economy, etc…) and what type of information you want to highlight to come up with a basis that next will be used to filter your signal and compress it.
If I find some time today I will try to make an example in Julia !
Forgot some references:
@wulpuqu , Thank you the response!!
I am using timeseries data like ECG and stocks. These kind of datasets can get messy or hard to work with. However, all the data is real valued and change in timely manner.
The goal is to change the shape of the signal to simpler form automatically, i.e. square to sine wave. I will use then this converted signal to perform my calculations.
Hence, the objective is to break down the signal (all the data) into set of simple smooth functions.
Now I can take individual filtered signal for calculations and later i can using statistics combine the individual results to obtain the final output.
{{{{{{{ The obtained output will be as if I worked on the original signal }}}}}}}
I will try to create a script based on your explanation. But please do send the example if you get time.
Thanks again!!
Here is an example. Taking the monthly mean temperature anomalies since 1880 as data.
I use the following script to filter the data:
using FFTW
using CSV
using Plots
## Some utils
# Wavelets
gmw(a,b,g) = begin
w_s = a*w
return normalize( (w_s.^b) .* exp.(-w_s.^g))
end
normalize(x) = x/sqrt(sum(abs2,x))
file_temp = "temp.csv"
if !isfile(file_temp)
run(`curl -L https://datahub.io/core/global-temp/r/1.csv -o $file_temp`)
end
d = reverse(CSV.File(file_temp).Mean) #Starting at year 1880
N = length(d)
d_fft = fft(d)
w = 2*pi*(0:N-1)/N
# low pass
a=18
b=0.1
g=1
psi_1=gmw(a,b,g)
N_zero = findfirst(y->y<1e-4,psi_1[2:end]) #starting at 2 since psi_1[1] is 0
conv_freq = psi_1 .* d_fft
conv_freq = vcat(conv_freq[1],conv_freq[2:N_zero],reverse(conj(conv_freq[2:N_zero])))
d_new = real(ifft(conv_freq))
N_css = (N_zero-1)*2+1
bad_N = div(N_css,2)+1
bad_N_css = round(Int,bad_N*N_css/N)
t = 1:N
t_new = range(1,N,length=N_css)
p_time=plot(t,d,label="d")
plot!(p_time,t_new[bad_N_css:N_css-bad_N_css],d_new[bad_N_css:N_css-bad_N_css],label="conv(d,psi)")
plot!(p_time,xaxis="Months",yaxis="°C",legend=:topleft)
p_fft=plot(w[1:div(N,2)+1],abs.(normalize(d_fft))[1:div(N,2)+1],label="fft(d)")
plot!(p_fft,w[1:div(N,2)+1],abs.(psi_1)[1:div(N,2)+1],label="psi")
plot!(p_fft,xaxis="Frequency (Pulsations)",yaxis="Amplitude",legend=:topleft)
p = plot(p_fft,p_time,layout=(2,1))
comp_factor = 100*(N-N_css)/N
println("Compression $comp_factor %")
# WARNING: working in Fourrier assume that all our signals are periodic, so we
# should pad left and right the signals in time when convoluting the signal with a filter. For simplicity, it is not done here !
# bad_N are the number of bad coefficients at each extremities
We get a compression of 84\%. Going from 3288 to 477 coefficients. Note that filtering signals induce a loss of information and some deformations. Here you can see that the filtered signal only coarsely match the trend of the initial data (not a correct padding here, see the next anwser of @rafael.guerra for properly filtering signals with DSP.jl).
With 100 Fourier Transform points, truncating the rest of the spectrum to zero and flippling input to avoid start/end discontinuity (no bells and whistles).
PS: edited the plot labels to better reflect the data displayed
using DSP, Plots, CSV, FFTW
file_temp = "temp.csv"
if !isfile(file_temp)
run(`curl -L https://datahub.io/core/global-temp/r/1.csv -o $file_temp`)
end
d = reverse(CSV.File(file_temp).Mean) #Starting at year 1880
d1 = [d; reverse(d)] # flip input to join start/end continuously
N = length(d1)
d1_fft = fft(d1)
t = 1:N;
f = LinRange(-0.5, 0.5, N+1)[1:N] # normalized frequency
fs = 1.0; # sampling frequency
fc = 0.01; # frequency cutoff, less than Fnyquist=0.5
filt = digitalfilter(Lowpass(fc,fs=fs), Butterworth(6));
d2 = filtfilt(filt, d1) #filtered signal
d2_fft = fft(d2)
Nc = 100 #Fourier spectrum samples kept (truncated beyond f >= Nc/N)
f0 = (Nc-1)/N # maximum frequency component kept
d3_fft_trunc = [d2_fft[1:Nc]; fill(0.0,N-2Nc+1); conj(d2_fft[Nc:-1:2])]
d3_trunc = real(ifft(d3_fft_trunc)) # time signal from truncated spectrum
p2 = plot(f, abs.(fftshift(d1_fft)), lw=2,lc=:blues,label="Flipped signal (2x3288 points)", xlim=(-0.1,0.1))
plot!(f, abs.(fftshift(d2_fft)), lw=2,lc=:reds,label="Filtered flipped signal",xlabel="Frequency")
plot!(f, abs.(fftshift(d3_fft_trunc)), ls=:dash,lw=1,lc=:black,label="Truncated spectrum (100 points)")
quiver!([-f0,f0],[180,180], quiver=([0,0],[-100,-100]), lc=:blue)
p1 = plot(t[1:NĂ·2], d1[1:NĂ·2], lc=:blues,label="Signal (3288 points)", legend=:topleft)
plot!(t[1:NĂ·2], d2[1:NĂ·2], lc=:reds, lw=3,label="Filtered flipped signal",xlabel="Time")
plot!(t[1:NĂ·2], d3_trunc[1:NĂ·2], ls=:dash,lc=:black, lw=1,label="From truncated spectrum (100 points)")
plot(p1,p2,layout=(2,1), size=(1000,1000))
@rafael.guerra Thank you for this wonderful example!!! Highly appreciate it.
@wulpuqu Thanks for the example and great explanation on the usage, highly appreciate it !!
