Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
How to find frequency

Subject: How to find frequency

From: Tomek

Date: 23 Mar, 2010 21:44:46

Message: 1 of 12

Hi,

I have a vector that has a periodic sequence of ones in it. For example
v = [0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 1 0 0 0 0 1 1 0 0 1]
I need to find the frequency of ones. However v may be disturbed either
by false ones that shouldn't be there or by ones that are slightly
shifted as in v(19). Also the number of ones in a group may vary. In v
there is usually a group of 2 ones followed by 5 zeros but number of
ones in a group my be different from 2. The run of zeros between groups
of ones (unless there is some false signal in between) is in my case in
the order of hundreds.

Thanks!

Subject: How to find frequency

From: Walter Roberson

Date: 23 Mar, 2010 22:04:43

Message: 2 of 12

Tomek wrote:

> I have a vector that has a periodic sequence of ones in it. For example
> v = [0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 1 0 0 0 0 1 1 0 0 1]
> I need to find the frequency of ones. However v may be disturbed either
> by false ones that shouldn't be there or by ones that are slightly
> shifted as in v(19). Also the number of ones in a group may vary. In v
> there is usually a group of 2 ones followed by 5 zeros but number of
> ones in a group my be different from 2. The run of zeros between groups
> of ones (unless there is some false signal in between) is in my case in
> the order of hundreds.

t = fft(v);
[unused, freq] = max(abs(real(t(2:floor(end/2)))));
usualrunlength = length(v) / freq;

In the case of your v, once you have dropped t(1) {which is the DC component},
the maximum absolute value in the first half of the fft occurs at position 4.
Your v is 28 elements long, and 28/4 = 7, which is exactly the length of your
typical run (0 0 0 0 0 1 1).

Subject: How to find frequency

From: us

Date: 23 Mar, 2010 22:11:22

Message: 3 of 12

Tomek <nomail@nonono.no> wrote in message <hobcof$283$1@news.onet.pl>...
> Hi,
>
> I have a vector that has a periodic sequence of ones in it. For example
> v = [0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 1 0 0 0 0 1 1 0 0 1]
> I need to find the frequency of ones. However v may be disturbed either
> by false ones that shouldn't be there or by ones that are slightly
> shifted as in v(19). Also the number of ones in a group may vary. In v
> there is usually a group of 2 ones followed by 5 zeros but number of
> ones in a group my be different from 2. The run of zeros between groups
> of ones (unless there is some false signal in between) is in my case in
> the order of hundreds.
>
> Thanks!

a hint:
- if i understand your problem correctly...

     help histc;
% eg,
     [nf,ix]=histc(v,0:1);

us

Subject: How to find frequency

From: TideMan

Date: 23 Mar, 2010 23:36:17

Message: 4 of 12

On Mar 24, 11:11 am, "us " <u...@neurol.unizh.ch> wrote:
> Tomek <nom...@nonono.no> wrote in message <hobcof$28...@news.onet.pl>...
> > Hi,
>
> > I have a vector that has a periodic sequence of ones in it. For example
> > v = [0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 1 0 0 0 0 1 1 0 0 1]
> > I need to find the frequency of ones. However v may be disturbed either
> > by false ones that shouldn't be there or by ones that are slightly
> > shifted as in v(19). Also the number of ones in a group may vary. In v
> > there is usually a group of 2 ones followed by 5 zeros but number of
> > ones in a group my be different from 2. The run of zeros between groups
> > of ones (unless there is some false signal in between) is in my case in
> > the order of hundreds.
>
> > Thanks!
>
> a hint:
> - if i understand your problem correctly...
>
>      help histc;
> % eg,
>      [nf,ix]=histc(v,0:1);
>
> us

Isn't that interesting?
Walter thinks by "frequency" the OP means the inverse of wave period,
whereas Urs thinks they mean the frequency of occurrence.
I think either could be right.
OP needs to clarify.

Subject: How to find frequency

From: Tomek

Date: 23 Mar, 2010 23:55:23

Message: 5 of 12

On 24.03.2010 00:36, TideMan wrote:
> On Mar 24, 11:11 am, "us "<u...@neurol.unizh.ch> wrote:
>> Tomek<nom...@nonono.no> wrote in message<hobcof$28...@news.onet.pl>...
>>> Hi,
>>
>>> I have a vector that has a periodic sequence of ones in it. For example
>>> v = [0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 1 0 0 0 0 1 1 0 0 1]
>>> I need to find the frequency of ones. However v may be disturbed either
>>> by false ones that shouldn't be there or by ones that are slightly
>>> shifted as in v(19). Also the number of ones in a group may vary. In v
>>> there is usually a group of 2 ones followed by 5 zeros but number of
>>> ones in a group my be different from 2. The run of zeros between groups
>>> of ones (unless there is some false signal in between) is in my case in
>>> the order of hundreds.
>>
>>> Thanks!
>>
>> a hint:
>> - if i understand your problem correctly...
>>
>> help histc;
>> % eg,
>> [nf,ix]=histc(v,0:1);
>>
>> us
>
> Isn't that interesting?
> Walter thinks by "frequency" the OP means the inverse of wave period,
> whereas Urs thinks they mean the frequency of occurrence.
> I think either could be right.
> OP needs to clarify.
>
Walter's approach is what I was looking for. I also need the
distribution of values but I do it simple way by sum(v) for ones and
length(v)-sum(v) for zeros. I think histc is too much for this simple case.

Subject: How to find frequency

From: Tomek

Date: 24 Mar, 2010 00:15:53

Message: 6 of 12

On 23.03.2010 23:04, Walter Roberson wrote:
> t = fft(v);
> [unused, freq] = max(abs(real(t(2:floor(end/2)))));
> usualrunlength = length(v) / freq;

Thanks Walter. I had similar idea but didn't take real() of fft. Why is
it necessary to do so?

Subject: How to find frequency

From: Tomek

Date: 24 Mar, 2010 01:00:02

Message: 7 of 12

On 23.03.2010 23:04, Walter Roberson wrote:
> Tomek wrote:
>
>> I have a vector that has a periodic sequence of ones in it. For example
>> v = [0 0 1 1 0 0 0 0 0 1 1 0 0 0 0 0 1 0 1 0 0 0 0 1 1 0 0 1]
>> I need to find the frequency of ones. However v may be disturbed
>> either by false ones that shouldn't be there or by ones that are
>> slightly shifted as in v(19). Also the number of ones in a group may
>> vary. In v there is usually a group of 2 ones followed by 5 zeros but
>> number of ones in a group my be different from 2. The run of zeros
>> between groups of ones (unless there is some false signal in between)
>> is in my case in the order of hundreds.
>
> t = fft(v);
> [unused, freq] = max(abs(real(t(2:floor(end/2)))));
> usualrunlength = length(v) / freq;
>
> In the case of your v, once you have dropped t(1) {which is the DC
> component}, the maximum absolute value in the first half of the fft
> occurs at position 4.
> Your v is 28 elements long, and 28/4 = 7, which is exactly the length of
> your typical run (0 0 0 0 0 1 1).

Hi again,

I've tried to run it on this sequence
TZ = 100;
v = [repmat([ones(1,7) zeros(1,130)],1,10) ones(1,7) zeros(1,TZ)];
and the computed runlength is extremely dependent on the number of
leading/trailing zeros (of course this number is shorter than typical
run length for zeros as in the example above). For TZ=100 I get
runlength of 34 (should be 137).
Is there any way to get the correct runlength even if the first/last
group of zeros is not equal to usual runlength for zeros?
The solution that comes to my mind is to remove all leading zeros as
well as trailing zeros and last group of ones. This way we'll only get
the part of v that is guaranteed to be roughly a multiple of
usualrunglength. But I'm not able to check whether the last group of
ones that I'm removing is not an accidental false peak until I find the
usualrunlength. If it would be a false peak then not all trailing zeros
would be removed and, again, usualrunlength would be incorrect.

Subject: How to find frequency

From: Walter Roberson

Date: 24 Mar, 2010 01:56:01

Message: 8 of 12

Tomek wrote:
> On 23.03.2010 23:04, Walter Roberson wrote:
>> t = fft(v);
>> [unused, freq] = max(abs(real(t(2:floor(end/2)))));
>> usualrunlength = length(v) / freq;
>
> Thanks Walter. I had similar idea but didn't take real() of fft. Why is
> it necessary to do so?

The imaginary part of the fft indicates phase information, which in
cases like these is more a reflection of how much "slippage" there is
around that frequency. On the other hand, when a peak is distorted by
phase, then it might be flattened in the real part, so perhaps just
taking abs() without the real() might work better.... Sorry, I'm still
fighting to understand the properties of fft myself.

I saw your test that you posted about later, and I will try it out, but
I'm getting some pretty slow reaction times in connecting from home to
work and displaying back to home, so real exploration might not be
practical until tomorrow.

Subject: How to find frequency

From: Tomek

Date: 24 Mar, 2010 02:35:10

Message: 9 of 12

On 24.03.2010 02:56, Walter Roberson wrote:
> Tomek wrote:
>> On 23.03.2010 23:04, Walter Roberson wrote:
>>> t = fft(v);
>>> [unused, freq] = max(abs(real(t(2:floor(end/2)))));
>>> usualrunlength = length(v) / freq;
>>
>> Thanks Walter. I had similar idea but didn't take real() of fft. Why
>> is it necessary to do so?
>
> The imaginary part of the fft indicates phase information, which in
> cases like these is more a reflection of how much "slippage" there is
> around that frequency. On the other hand, when a peak is distorted by
> phase, then it might be flattened in the real part, so perhaps just
> taking abs() without the real() might work better.... Sorry, I'm still
> fighting to understand the properties of fft myself.
Thank you very much. I realized what I did wrong in my initial trial. I
forgot to do the last division of length(v) / freq. Instead I thought
freq was the solution. As far as I see, the abs()-only version would
also work. I'm not sure which one is more appropriate.

>
> I saw your test that you posted about later, and I will try it out, but
> I'm getting some pretty slow reaction times in connecting from home to
> work and displaying back to home, so real exploration might not be
> practical until tomorrow.
I think the reason might be the resolution of fft. In the signal
spectrum there are very strong peaks that are not reflected in fft
spectrum because these strong peaks occur somewhere between points were
frequency spectrum is sampled. However increasing resolution of fft is
not a simple fix for this. Once again, thank you for your help.

Subject: How to find frequency

From: TideMan

Date: 24 Mar, 2010 02:53:02

Message: 10 of 12

On Mar 24, 3:35 pm, Tomek <nom...@nonono.no> wrote:
> On 24.03.2010 02:56, Walter Roberson wrote:> Tomek wrote:
> >> On 23.03.2010 23:04, Walter Roberson wrote:
> >>> t = fft(v);
> >>> [unused, freq] = max(abs(real(t(2:floor(end/2)))));
> >>> usualrunlength = length(v) / freq;
>
> >> Thanks Walter. I had similar idea but didn't take real() of fft. Why
> >> is it necessary to do so?
>
> > The imaginary part of the fft indicates phase information, which in
> > cases like these is more a reflection of how much "slippage" there is
> > around that frequency. On the other hand, when a peak is distorted by
> > phase, then it might be flattened in the real part, so perhaps just
> > taking abs() without the real() might work better.... Sorry, I'm still
> > fighting to understand the properties of fft myself.
>
> Thank you very much. I realized what I did wrong in my initial trial. I
> forgot to do the last division of length(v) / freq. Instead I thought
> freq was the solution. As far as I see, the abs()-only version would
> also work. I'm not sure which one is more appropriate.
>
>
>
> > I saw your test that you posted about later, and I will try it out, but
> > I'm getting some pretty slow reaction times in connecting from home to
> > work and displaying back to home, so real exploration might not be
> > practical until tomorrow.
>
> I think the reason might be the resolution of fft. In the signal
> spectrum there are very strong peaks that are not reflected in fft
> spectrum because these strong peaks occur somewhere between points were
> frequency spectrum is sampled. However increasing resolution of fft is
> not a simple fix for this. Once again, thank you for your help.

IMHO, using FFT is the wrong approach.
You should be looking for jumps:
dv=diff(v);
Now, in dv, the +1s are an upward jump and the -1s are a downward
jump.
up=find(dv==1);
down=find(dv==-1);
Now you can play with the addresses up and down to get the interval
between spikes.
The width of each spike is: down - up
The interval between spikes is: diff(up)
etc
etc

Subject: How to find frequency

From: Walter Roberson

Date: 24 Mar, 2010 05:26:17

Message: 11 of 12

Tomek wrote:

> I've tried to run it on this sequence
> TZ = 100;
> v = [repmat([ones(1,7) zeros(1,130)],1,10) ones(1,7) zeros(1,TZ)];
> and the computed runlength is extremely dependent on the number of
> leading/trailing zeros

You are right; I have done some testing, and finding the peak fft
coefficient is *not* a good way of detecting the period of a binary
signal. For example, [0 ones(1,63)] is completely flat (and real-valued)
in the fft after the DC component, all -1's, which would mean
-sin(x)-sin(2*x)-sin(3*x)-sin(4*x) and so on. [ones(1,63) 0] on the
other hand is similar to a sine curve, peaking in the middle, starting
from negative coefficients and going to positive coefficients. As
humans, we can recognize that how-ever we want to define the period of
such a vector to be, then by symmetry it should be the same in the two
cases.

Each zero padding on the end of data adds one more sine convolution (of
the fft coefficients). And repeating the signal, fft([x x]) can have the
effect of moving the detected peak from location n to location 2*n+1...
Oh wait, that one actually makes sense since you doubled the length so
length / location would still be either the same or a hair more than
before (increased accuracy.)

Anyhow, I do not have any fft-based suggestions on how to proceed.
Tideman hinted at other methods; going along with what he was saying,
one approach might be to histogram the detected periods. If you get a
notable multi-modal distribution in the histogram, that would be a
reflection of internal transitions within the periodic binary signal...
e.g., [1 1 0 1 0 0 0 0 0] repeated (without error) many times would have
as many bins with "gap length 1" as for "gap length 5". Might be a bit
tricky to reconstruct the periodic signal from the bin counts... for
example, 0 1 0 1 0 0 0 0 1 would have twice as many zeros-gaps of length
1 as it would of length 4. I guess you could look for the longest gap
runs of zeros that have significant counts, and then look for other bins
that have close approximations of multiples of that count.. but that
won't tell you the order, and you will need to adapt this (perhaps by a
second histogram counting runs of 1's) to figure out the likely total
signal length.


This signal... when there are varying runs of 1's, is there a known (or
reasonably estimated) probably distribution for the number of 1's?
Poisson perhaps? Or does it represent something like "work load varying
over time in an way whose prediction is beyond the scope of this
project"? Or is it just random measurement error, or relatively accurate
measurement of a process that is prone to minor errors?

Subject: How to find frequency

From: David Young

Date: 24 Mar, 2010 08:41:06

Message: 12 of 12

Walter Roberson <roberson@hushmail.com> wrote in message <hobrfi$nco$1@canopus.cc.umanitoba.ca>...
> ...
>
> The imaginary part of the fft indicates phase information, which in
> cases like these is more a reflection of how much "slippage" there is
> around that frequency. On the other hand, when a peak is distorted by
> phase, then it might be flattened in the real part, so perhaps just
> taking abs() without the real() might work better....
> ...

Just abs() is much better.

The imaginary and real parts of the fft contribute equally to the phase, and also contribute equally to the magnitude. It's not the case that the imaginary part carries any more phase information than the real part: phase is atan2(imagpart, realpart).

Thus you could have a peak in the spectrum that happened to have a phase of pi/2, but the real part of the fft would be zero at that frequency - the peak would only show up in the imaginary part, and you'd miss it if you used real(). However, abs() will show a peak regardless of the phase.

So, if you want to ignore the phase, and look at the magnitude or energy of the signal, just use abs.

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us