Hello,
I am looking for a way to determine the instantaneous frequency of a signal (filtered BDLC motor phase voltage) to determine the speed of this motor.
The sampling frequency is 8 kHz, the max, electrical frequency about 1 kHz.
For Python I found the following code: https://gist.github.com/endolith/255291
Any idea, how to do this in Julia? There is quite some scientific literature on this topic, but I did not found any algorithm that is straightforward to implement. I would like to reach a time resolution of at least 10 ms and an accuracy of at least 1% for the frequency. This is not easy, taking the low sampling frequency into account.
On your original links, you could try to convert those to Julia. Most don’t look that bad. I’d use the following version rather than the Gist, because it has an explicit license (MIT).
I know this thread has been inactive, but I’ve been learning Julia and having read this thread I thought I would try to convert a frequency estimator from python.
So, here’s a simple algorithm that does pretty well at finding the frequencies.
Note, it forms a big matrix from the input array and takes an SVD, so, maybe start with a subset of the vector first.
"""
esprit(x::Array, M::Int, p::Int, Fs::Float64=1.0)
ESPRIT algorithm for frequency estimation.
Estimation of Signal Parameters via Rotational Invariance Techniques
Given length N signal "x" that is the sum of p sinusoids of unknown frequencies,
estimate and return an array of the p frequencies.
# Arguments
- `x::Array`: complex length N signal array
- `M::Int`: size of correlation matrix, must be <= N.
The signal subspace is computed from the SVD of an M x (N-M+1) signal matrix
formed from N-M+1 length-M shifts of the signal x in its columns.
For best performance for 1 sinusoid, use M = (N+1)/3 (according to van der Veen and Leus)
For faster execution (due to smaller SVD), use small M or small N-M
- `p::Int`: number of sinusoids to estimate.
- `Fs::Float64`: sampling frequency, in Hz.
# Returns
length p real array of frequencies in units of Hz.
"""
function esprit(x::Array, M::Int, p::Int, Fs::Float64=1.0)
N = length(x)
X = x[ (1:M) .+ (0:N-M)' ]; # in julia, likely faster to loop than to allocate huge index array
U,s,V = svd(X);
D,_ = eig( U[1:end-1,1:p] \ U[2:end,1:p] );
angle.(D)*Fs/(2*π)
end
I’m steadily converting a lot of little utility Matlab functions to Julia. For a complex time series, this function should give a good estimate of the instantaneous frequency:
"""
USAGE: f_Hz = instfreq(t_s,z)
"""
function instfreq(t,z)
ph = unwr(angle(z));
dt = t[2]-t[1];
omeg = cdif(ph)/dt;
f = omeg/2/pi;
return f;
end
along with these two other functions:
"""
USAGE: pfix = unwr(praw,base{2*pi} [,ref])
"""
function unwr(praw,base=2π,ref=nothing)
m = length(praw);
pfix = praw;
if(ref ≢ nothing) # rewrap to reference phase(s) provided
pfix .-= base*round((pfix-ref)/base);
return pfix;
end
p = praw;
pd = diff(p);
ijump = find(abs.(pd) .> base/2);
for i in 1:length(ijump)
k = ijump[i];
p[(k+1):m] .-= base * round(pd[k]/base);
end
pfix = p;
return pfix;
end
and
"""
USAGE: ydot = cdif(y)
"""
function cdif(y)
n = length(y);
ydot = [diff(y[[1,3]]);
y[3:n]-y[1:(n-2)];
diff(y[[n-2,n]])]/2;
return ydot;
end