MATLAB Answers

## Detection of sinusoids in the discrete sequence

Asked by Andrey

### Andrey (view profile)

on 5 Oct 2012

Good day! Please excuse me for my English. I have the following problem: There is a period of a discrete signal. This signal is the sum of sinusoidal signals (sinusoids). These sinusoids not always begin and end at the beginning and end of the period. As an example, I suggest the following (see figure) Period length = 300. There are two sinusoidal signals: one signal all over the period (300), the second - only half of the period (150). Question: how to detect both signals in this example? That is, how to determine the frequency of these two signals. If both signals were all over the period, the usual Fourier (fft) would perform this task. But in this example, Fourier did not determine the frequency correctly. Thank you in advance! http://i064.radikal.ru/1210/90/ed1741e392a6.jpg

Greg Heath

### Greg Heath (view profile)

on 5 Oct 2012

The link doesn't work. I get the error message

The page you were looking for doesn't exist.

You may have mistyped the address or the page may have moved.

Will you please submit a short piece of code that will generate the signal and produce the plot?

TIA

Andrey

### Andrey (view profile)

on 8 Oct 2012

Sorry. I corrected

## Products

No products are associated with this question.

## 4 Answers

### Greg Heath (view profile)

Answer by Greg Heath

### Greg Heath (view profile)

on 6 Oct 2012

If both sampling frequencies are equal to Fs = 1/dt, then the folowing conditions should be met

1. The signal amplitudes are comparable so that neither spectrum peak is comparable to the other signal spectrum's peak sidelobe.

2. Each signal length is longer than a fundamental period:

t1max > 1/f1 and t2max = 2*t1max > 1/f2

3. If x1 is zeropadded to length t2max, the peak spectra can be resolved when

df = 1/(t2max+dt) < abs(f1-f2)

Thank you for formally accepting my answer.

Answer by Andrey

on 8 Oct 2012
Edited by Andrey

### Andrey (view profile)

on 8 Oct 2012

Thanks for the answers! This example was very simple. In this example I would draw your attention to the following: If the sinusoid is not defined on the whole window and when it intersects with other unknown sinusoid, the Fourier is not exactly determine the frequency. Of course, we can take a smaller window, but .... How about this example: Note that I am defining period of sinusoid (not the frequency). That is, how many points should be sinusoid to make the full period. Sorry, but I feel so comfortable. In my example, sinusoids with three different lengths, with different periods (frequencies). It are put together and create a new signal. Objective: to determine the period (frequency) of the source sinusoids.

function [ ] = Test()
x = 1:100;
y1 = CreateSin(x, 20, 78, 30);
y2 = CreateSin(x, 13, 60, 45);
y3 = CreateSin(x, 3, 100, 60);
y = y1 + y2 + y3;
subplot(4, 1, 1);
plot(x, y1);
axis([0, 100, -2, 2]);
subplot(4, 1, 2);
plot(x, y2);
axis([0, 100, -2, 2]);
subplot(4, 1, 3);
plot(x, y3);
axis([0, 100, -2, 2]);
subplot(4, 1, 4);
plot(x, y);
axis([0, 100, -2, 2]);
end
% function create sinusoid
% x - domain
% begin - point from which to begin a sinusoid
% finish - the point at which ends sinusoid
% period - sinusoid period
function [ y ] = CreateSin(x, begin, finish, period )
y = zeros(1, length(x));
y(begin:finish) = y(begin:finish) + sin(x(begin:finish) * 2 * pi / period);
end

This example is not so simple. Because there is sinusoids is almost completely overlap. And sinusoids very short. This makes the task very difficult. I would be grateful for any ideas

Answer by Andrey

on 24 Oct 2012
Edited by Andrey

### Andrey (view profile)

on 24 Oct 2012

Here is a fantastic example:

Imagine a planet. The water level at the beach depends on the moon. We all know that the moon pulls the water. And on the last day of our planet 100 Earth hours. Here's the plot: http://s019.radikal.ru/i611/1210/2a/7ba1b19ff15d.jpg

Now imagine that we have a second moon, which rotates twice as fast around the planet (the period is 50 Earth hours) Moons can be at one point, or can on the different sides of the moon. Therefore, the gravitational force of the two moons are summarized as two vectors, so we get the following picture. http://s005.radikal.ru/i210/1210/48/37f102edf9bd.jpg

And that's not all. Second Moon permanently disappear somewhere. No one knows where. Just disappears. So when the second moon is absent the effect of gravity on the water only to the first moon. Our second moon appears stable at 27 o'clock, and disappears by 64 o'clock http://s14.radikal.ru/i187/1210/d2/86c0e6daae87.jpg

Add to our skies third moon: http://s018.radikal.ru/i500/1210/5c/158fc57c9994.jpg

Now imagine that we want to know about when the moon disappear. We also need to know how fast to spin the moon around the planet. And all this we can know, having the only the graph (the water level):

Who has ideas how to do this algorithmically?

function [ ] = Run()
close all;
clear all;
clc;
moon1 = CreateMoon(100, 1, 100, 100);
moon2 = CreateMoon(100, 27, 64, 50);
moon3 = CreateMoon(100, 50, 70, 75);
subplot(4, 1, 1);
plot(moon1, '.');
axis([0, 100, -2, 2]);
subplot(4, 1, 2);
plot(moon2, '.');
axis([0, 100, -2, 2]);
subplot(4, 1, 3);
plot(moon3, '.');
axis([0, 100, -2, 2]);
subplot(4, 1, 4);
plot(moon1+moon2+moon3, '.');
axis([0, 100, -2, 2]);
end
% function create moon
% dayLength - the number of hours in a day
% appearance - time of appearance
% disappearance - time of disappearance
% period - time for a complete orbit around the planet
function [ y ] = CreateMoon(dayLength, appearance, disappearance, period )
x = 1:dayLength
y = zeros(1, length(x));
y(appearance:disappearance) = y(appearance:disappearance) + sin(x(appearance:disappearance) * 2 * pi / period);
end

Answer by Andrey

on 24 Oct 2012
Edited by Andrey

### Andrey (view profile)

on 25 Oct 2012

I am interested in a method that, in my view, better helps to solve my problem: Empirical Mode Decomposition (EMD) Look what results achieved by the example with moons http://s56.radikal.ru/i154/1210/cd/322b7fd4904b.jpg

Unfortunately, we were able to identify only one moon (which is the full width of the window). But we were able to algorithmically find at least one of the signal composes. I think this is a good start to finding a solution.

moon1 = CreateMoon(1000, 1, 1000, 100);
moon2 = CreateMoon(1000, 340, 1000, 33);
moon3 = CreateMoon(1000, 1, 700, 70);

http://s017.radikal.ru/i442/1210/f8/235bf05cd0a4.jpg

Look at the results and compare them to the sine wave. Of course, there is no 100% accurate. But we're just getting started. I marked in red pencil, which would have placed the boundary of the sine wave and their periods.

Can yCan you tell me something about EMD (or other similar methods) that would help me move on?ou tell us anything about the EMD, that would help me move on?

#### Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply today