See http://gist.github.com/255291 for more methods
-
-
Save endolith/129445 to your computer and use it in GitHub Desktop.
function [frequency, period, crossings, seconds, samples]=freq(wavefile, siz) | |
% FREQ Measures frequency of mono .WAV files by counting zero crossings | |
% | |
% [frequency, period, crossings, seconds, samples]=freq(wavefile, siz, fmt24bit) | |
% | |
% frequency : in hertz | |
% period : in seconds | |
% crossings : the total number of zero crossings | |
% seconds : total amount of time between first and last zero crossing | |
% samples : total number of samples between first and last zero crossing | |
% | |
% wavfile : file name | |
% (The ".wav" extension is appended if no extension is given.) | |
% siz : 'size' returns the size of the audio data | |
% : [N1 N2] returns only samples N1 through N2 | |
% : N1 returns the first N1 samples | |
% : default is all the data | |
% endolith@gmail.com July 12th, 2005 | |
% Meant to measure frequency very accurately, for comparing ADC clocks. Only useful in the specific circumstance I was using it, with long wav files and signal much larger than noise, so that noise doesn't cause spurious zero crossings | |
% Show help | |
if nargin < 1 | |
help freq; | |
return; | |
end | |
%Load WAV file | |
if nargin < 2 | |
[sig, sf, bits] = wavread (wavefile); | |
else | |
if nargin < 3 | |
[sig, sf, bits] = wavread (wavefile, siz); | |
end | |
end | |
%We need high precision for these numbers | |
output_precision(9) | |
%Find all zero crossings | |
c=[(sig(1:size(sig,1)-1)<=0).*(sig(2:size(sig,1))>0);0]; | |
%Find locations of the rising edge zero crossings | |
a=find(c); | |
%Calculate number of samples between first and last zero crossing (relevant data) | |
samples=(a(size(a,1)))-a(1) | |
% Should specify the precision of the measurement (1000.0 Hz for a short sample, 1000.00000 Hz for a long sample?) | |
%Calculate time in seconds of relevant section | |
seconds=samples/sf | |
%Number of crossings in relevant section | |
crossings=size(a,1)-1 | |
%Frequency is crossings per second | |
frequency=crossings*sf/samples | |
period=1/frequency | |
%Plot | |
%closeplot | |
%d=zeros(size(sig,1),1); | |
%d(a(1):a(size(a,1))-1)=1; | |
%dsize = size(d(d>0),1) | |
%t = [1:size(sig,1)]; | |
%plot(t,sig,t,c,t,d) | |
%allcs = size(c(c>0),1) | |
from __future__ import division | |
from numpy import logical_and, average, diff | |
from matplotlib.mlab import find | |
def freq_from_crossings(sig, fs): | |
"""Estimate frequency by counting zero crossings | |
Doesn't work if there are multiple zero crossings per cycle. | |
""" | |
indices = find(logical_and(sig[1:] >= 0, sig[:-1] < 0)) | |
# Linear interpolation to find truer zero crossings | |
crossings = [i - sig[i] / (sig[i+1] - sig[i]) for i in indices] | |
return fs / average(diff(crossings)) |
@LaurentNorbert Yeah, that's right, though it's only going to take microseconds either way.
average(diff(crossings))
can be replaced by
(crossings[-1]-crossings[0])/(len(crossings)-1)
The average method would allow you to do a trimmed mean to throw away outliers, etc.
I'm new to python. Can someone tell me what sig and fs mean? Thanks.
@lewham97 those are just variable names, not specific to python. sig = signal, fs = sample rate
@endolith Thanks for clearing that up. Would this code be suitable for calculating voltage zero crossings from the output of a comparator circuit?
@lewham97 I'm not sure what you mean. It counts zero crossings of a constant frequency and calculates the frequency
@endolith I’m looking to count zero crossings of a square wave with a peak of 3.3v using a Raspberry Pi to calculate the frequency. I’m hoping that this code will be suitable for the job!
I'm amazed by the latter Python code which doesn't seem optimized at all.
Don't we have sum(diff(crossings)) = crossings[-1] - crossings[0]?
By the way, I used this function on an Arduino, in C. Of course, the reference point is not 0. but 512 if the signal is quantized on 10 bits.