Count zero crossings to determine frequency

readme.md

See http://gist.github.com/255291 for more methods

freq.m

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


% [email protected] 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)

freq_from_crossings.py

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))

Indices = find(...) is always empty.

import numpy as np
indices = np.nonzero(np.diff(sig > 0))[0]

This might work!

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.

@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!