Detection of sinusoids in the discrete sequence

[included]

A chirp is a small burst of a signal somewhere inside the sample period, it could be a burst of a single frequency or an increasing frequency signal or even an impulse.

When you define your analysis block B2 over the sampled block B1 you only need to consider the data in B1 that coincides with B2 and all the data outside is not included in the analysis.

Say you have a sample block 1000 points long, you want to analyse the data from point 200 through to point 299 so just that part is extracted and processed and the rest of the data is ignored.

Once you have defined the shape of your chirp you can then design quicker ways to do tha analysis. If it is a single frquency then a digital filter can be applied and as you have desrcibed your problem it appears you are indeed looking for single frequency chirps.

You could use a faster way to hone in on your key frequencies and then apply the rigorous method to just those specific frequencies, that should take less than an hour - maybe tens of minutes. Once you find the good solution you probably should code it in macro to get it to run fast enough.

 

[S-LifeLover]

Now overlay B2 on to of B1 so the B2 starts at the first point in B1 and multiply the blocks and sum the result. Store it. Just sum the result of B1*B2 and exclude the rest of B1.

B1*B2.
^
You mean the elementwise multiplication or finding correlation of two signals B1 and B2?
Thanks.

 

[included]

---

You can do a correlation but the only way to extract the signals you want is to do it the hard way, elementwise and summate.

The faster way is to use your original 'query signal' ie B2 and do a fourier domain correlation but instead calculate the coherance, it is fast and will give you the best choice of your 'query signal', then go back to the hard way to extract as much data as you can.

But now you dealing with 3 blocks of data so you need to keep track of what you do. I would start with an FFT base line and get a coherance map for your signal, then go to DFT base to increase accuracy once you have a peak in the coherance. Single signal coherant behaviour was used 30 years ago to detect non linear damping, rattles, frequency migration etc and so is not new but better you get your core code for the block manipulation started.

The DFT method could replace basic block manipulation but speed issues are unknown as it all depends on your sample size. I do not even know how big your sample sequence is and its real time rate, (sample frequency)

The reason I answer your post is I did this for many years but if I write the code for you - who learns to interpret the results. If you do not do it yourself it is worthless.

 

[S-LifeLover]

Good day!
I tried to write the code of your algorithm. Could you take a look at this code?

clear all;
close all;
clc;

length = 100;
period = 20;
y = CreateMoon(length, 1, 100, period);

Pereborka(y);

function [ y ] = sinus( len, phase, amp, per)
    radiansPerDot = ( pi * 2 ) / per;
    phase = phase * radiansPerDot;
    radians = (0:len-1) * radiansPerDot;
    y = amp * sin(radians + phase);
end

function [ ] = Pereborka( y )

    B1 = y;
    
    max_MultySum = 0;
    max_Pos = 0;
    max_Phase = 0;
    max_Len = 0;
    
    wb = waitbar(0, 'Working...');
    
    valuesLength = length(y);
    % curPosB2 - current position of B2
    % curPerB2 - current period of B2
    % curAmplB2 - current amplitude of B2
    % curPhaseB2 - current phase of B2
    curAmplB2 = 1;
    
    %P4
    for curPerB2 = 1 : valuesLength
        waitbar(curPerB2 / valuesLength, wb);
        
        %P3
        for curLenB2 = 1 : valuesLength
            %P2
            for curPhaseB2 = 1 : valuesLength
                %P1
                for curPosB2 = 1 : valuesLength
                    if(curPosB2 + curLenB2 - 1 > valuesLength)
                        break;
                    end;
                    curB2 = sinus( curLenB2, curPhaseB2, curAmplB2, curPerB2);
                    part_B1 = B1(curPosB2:curPosB2+curLenB2-1);
                    multySum = sum(curB2.*part_B1);
                    if(max_MultySum < multySum)
                        max_MultySum = multySum;
                        max_Pos = curPosB2;
                        max_Phase = curPhaseB2;
                        max_Len = curLenB2;
                    end;
                end;
            end;
        end;
    end;
    
    delete(wb);
    
    disp(['max_MultySum = ', num2str(max_MultySum)]);
    disp(['max_Pos = ', num2str(max_Pos)]);
    disp(['max_Phase = ', num2str(max_Phase)]);
    disp(['max_Len = ', num2str(max_Len)]);

    B2 = sinus( max_Len, max_Phase, curAmplB2, curPerB2);
    figure()
    hold on;
    plot(B1);
    plot(max_Pos:max_Pos+max_Len-1, B2, 'Color', 'Red');
    title(num2str(max_MultySum));

end

 

[included]

As I said earlier I do not code in Matlab so I will have to get a copy and read the manual to understand it fully. But first glance is a set of nested loops to change position, freq, amp and phase, then store the results in an array. Seems correct procedure but exact code ??? I will have to check.

Can you plot your B1 so we can see the original signal and then plot your output arrays. That would be usefull.

If you are using a sine chirp then this is equivent to taking the sin*cos product of B1 in the B2 overlay, this is faster in some instances but needs more interpretation and then further analysis, the hard way extracts all the data in one pass.

Of course if your chirp has a varying frequency sinusoid and/or damping the only way is the hard way.

Please post the results.

- - - Updated - - -

OK I have a manual and the process looks good, I dont beleive that changing amplitude is neceassry, just use a typical level as all it will do is scale the results.
Next I would put the results into a matrix Result(position,Length) and you could use two matices for the frequency and the phase or combine it as a complex number. The complex number could be further manilulated by testing the phase and then calculating the magnitude and just display that.
Then plot a SURF of the matrix - you should see one peak somewhere on the 3D plot.

But you still do not know how good a guess your chirp was so this is where you need to calculate the coherance between the guess chirp (B2) and the B1 at the location, length, Freq, Phase that gives you the peak reading.

You have all the information available to calculate Cxy from abs[CSxy]^2/B1xx*B2yy and for fun you could run this particular B2 across the whole of B1 again and as it is just one pass it should be quick.

Let me know how that goes.

 

[S-LifeLover]

Result:


As we can see, the result is not what we expected. The algorithm does not find a sine wave with a period of = 20. The algorithm found a sinusoid with period = 100.
I think we should use instead of the correlation sum. I'm working on it.


Blue line - source sinusoid.
Red line - result sinusoid.

OK I have a manual and the process looks good, I dont beleive that changing amplitude is neceassry, just use a typical level as all it will do is scale the results.

I did not change the amplitude, because I could not wait to end the program. But to solve my problem, I must calculate the amplitude.

 

[included]

The method does work and so somewhere in the code there is a glitch but as I am not well practised in Matlab I cant see it yet. My version works fine.

You should store all the multysum results into an array and then look at the array and see where things go wrong

The reason to do this way is you need to get the exact phase of B2 correct. Correlation is fine if use lots of B2 phase shifted blocks but then the result is the same thing as you have already done. But you need the core code of the slider so you can add the coherance check later.

 

[S-LifeLover]

You should store all the multysum results into an array and then look at the array and see where things go wrong

If period = 100 then
sum(B1.*B2) = 50

if period = 20 then
sum(B1.*B2) = 50

In fact, both versions give the same result. But on the chart there is only one. Therefore, such a result is obtained.
I'm working on keeping all sinusoids in the matrix. And to build the surfaces.

 

[included]

Result:
I did not change the amplitude, because I could not wait to end the program. But to solve my problem, I must calculate the amplitude.

Once you find the block of data where the signal is maximal to your chirp then simply do a new calculation on that specific data set to get the amplitude.

- - - Updated - - -

If period = 100 then
sum(B1.*B2) = 50

if period = 20 then
sum(B1.*B2) = 50

In fact, both versions give the same result. But on the chart there is only one. Therefore, such a result is obtained.
I'm working on keeping all sinusoids in the matrix. And to build the surfaces.

This cannot be, a rectified sine wave does not have an average value of 0.5 per bin, ie 50/100=0.5

A triangular wave maybe but not a sine wave.

Check the command is correct in Matlab

 

[S-LifeLover]

I apologize for the doubts in your algorithm. I fixed a couple of bugs in the program. Now the result is correct.
I continue to work on the program.

function [ ] = Pereborka( y )

    B1 = y;
    
    max_MultySum = 0;
    max_Pos = 0;
    max_Phase = 0;
    max_Len = 0;
    max_Per = 0;
    
    wb = waitbar(0, 'Working...');
    
    valuesLength = length(y);
    % curPosB2 - current position of B2
    % curPerB2 - current period of B2
    % curAmplB2 - current amplitude of B2
    % curPhaseB2 - current phase of B2
    curAmplB2 = 1;
    
    %P4
    for curPerB2 = 1 : valuesLength
        waitbar(curPerB2 / valuesLength, wb);
        
        %P3
        for curLenB2 = 1 : valuesLength
            %P2
            for curPhaseB2 = 1 : valuesLength
                %P1
                for curPosB2 = 1 : valuesLength
                    if(curPosB2 + curLenB2 - 1 > valuesLength)
                        break;
                    end;                    
                    
                    curB2 = sinus( curLenB2, curPhaseB2, curAmplB2, curPerB2);
                    part_B1 = B1(curPosB2:curPosB2+curLenB2-1);
                    multySum = sum(abs(curB2.*part_B1));
                    if(max_MultySum < multySum || (max_MultySum == multySum && max_Len < curLenB2))
                        max_MultySum = multySum;
                        max_Pos = curPosB2;
                        max_Phase = curPhaseB2;
                        max_Len = curLenB2;
                        max_Per = curPerB2;
                    end;
                end;
            end;
        end;
    end;
    
    delete(wb);
    
    disp(['max_MultySum = ', num2str(max_MultySum)]);
    disp(['max_Pos = ', num2str(max_Pos)]);
    disp(['max_Phase = ', num2str(max_Phase)]);
    disp(['max_Len = ', num2str(max_Len)]);

    B2 = sinus( max_Len, max_Phase, curAmplB2, max_Per);
    figure()
    hold on;
    plot(B1);
    plot(max_Pos:max_Pos+max_Len-1, B2, 'Color', 'Red');
    title(num2str(max_MultySum));

end

 

[included]

I ran the 20 period and the 100 period cases in sigproc and found one sums to 69 whilst the other sums to 0.

The plot below is the result after B1*B2 and is as one would expect. One is a fully rectified sine whilst the other mirrors about the x axis.

Glad the code is now finding this.

 

[S-LifeLover]

slightly modified athe code

function [ ] = bruteforce( y )

    B1 = y;
    
    max_MultySum = 0;
    max_Pos = 0;
    max_Phase = 0;
    max_Len = 0;
    max_Per = 0;
    
    wb = waitbar(0, 'Working...');
    
    valuesLength = length(y);
    % curPosB2 - current position of B2
    % curPerB2 - current period of B2
    % curAmplB2 - current amplitude of B2
    % curPhaseB2 - current phase of B2
    curAmplB2 = 1;
    
    %P4
    for curPerB2 = 1 : valuesLength
        waitbar(curPerB2 / valuesLength, wb);
        
        %P3
        for curLenB2 = 1 : valuesLength
            %P2
            for curPhaseB2 = 1 : valuesLength
                %P1
                for curPosB2 = 1 : valuesLength
                    if(curPosB2 + curLenB2 - 1 > valuesLength)
                        break;
                    end;                    
                    
                    curB2 = sinus( curLenB2, curPhaseB2, curAmplB2, curPerB2);
                    part_B1 = B1(curPosB2:curPosB2+curLenB2-1);
                    multySum = sum(curB2.*part_B1);
                    if(max_MultySum < multySum || (max_MultySum == multySum && max_Len < curLenB2))
                        max_MultySum = multySum;
                        max_Pos = curPosB2;
                        max_Phase = curPhaseB2;
                        max_Len = curLenB2;
                        max_Per = curPerB2;
                    end;
                end;
            end;
        end;
    end;
    
    delete(wb);
    
    disp(['max_MultySum = ', num2str(max_MultySum)]);
    disp(['max_Pos = ', num2str(max_Pos)]);
    disp(['max_Phase = ', num2str(max_Phase)]);
    disp(['max_Len = ', num2str(max_Len)]);
    disp(['max_Per = ', num2str(max_Per)]);

    B2 = sinus( max_Len, max_Phase, curAmplB2, max_Per);
    figure()
    hold on;
    plot(B1);
    plot(max_Pos:max_Pos+max_Len-1, B2, 'Color', 'Red');
    title(num2str(max_MultySum));

end

Good day!
Your algorithm is very good at finding a sine wave signal. Unfortunately, for my task this algorithm is not suitable.

Here is my example:

clear all;
close all;
clc;

y = CreateMoon(70, 1, 100, 35);
y = CreateMoon(70, 1, 100, 10);

bruteforce(y);

And result:

max_MultySum = 36.5255
max_Pos = 1
max_Phase = 4
max_Len = 70
max_Per = 38

Alas ... The algorithm was not able to help me in the case where the source signal consists of several sinusoids.
Perhaps I understand the algorithm not true.

 

[included]

Now subtract the signal from moon and do the process again to find the second frequency.

It looks pretty close to me but all depends on the resolution of moon. The method is very exact but does require some manipulation and code tweaks to make it work well.

What I do not understand is my code works very very well but yours seems to be of a little bit. Cant understand why but you calculate one maxsum whilst I store the whole array and look at the peaks so I see all the peaks and you only get one.

Try to store an array of results till you sort it out.

 

[S-LifeLover]

Now subtract the signal from moon and do the process again to find the second frequency.

Why? The program has found incorrect sinusoid. And if I subtract it from the original signal and run the program again, I get a second incorrect sinusoid.

The method is very exact but does require some manipulation and code tweaks to make it work well

Could you tell me about these tweaks?

Cant understand why but you calculate one maxsum whilst I store the whole array and look at the peaks so I see all the peaks and you only get one

I think that we need the maximum peak. And it usually a single.
And I'm not yet poor understanding how to make a visual representation of the array. Each sum depends on four parameters: the position, length, phase and period. And to build surface we need 2 parameters: x and y.
You suggested the following:

Next I would put the results into a matrix Result(position,Length) and you could use two matices for the frequency and the phase or combine it as a complex number.

But how do I calculate the phase of the length and position?
I count sum as follows:
sum = sum(B1.* BuildSinus(Length, Position, Phase, Period) )
That is, for each position, and the length, I have several phases.
In other words in my system equation 5 variables. And to build the plane in the system of equations must be 3 variables.

 

[included]

I do not know why your code calculates the wrong frequency, something in there is not quite right that is why I suggest to store is an array.

A(x,y,z) where x is the position, y is the length and z is the complex representation of B1 for those parameters.
You could if you wish make A(x,y,z) have z = period and A1(x,y,w) have w = phase.

Now Surf A(...) and A1 and see what is going on. The array A will of course be skewed as less points are stored as the length is increased. This suface describes your analysis completely.

Tweaks. Well in development I would have a block size of 128 and also do a FFT on it to see what that looks like and ascertain basic signal levels for the peaks. Those would be used in the slider you amusingly called bruteforce.
Or a DFT for odd block sizes.

Then I would also do a coherance check at the frequencies that come out of the code to see how good the fit is. Using a spectral map is the easiest but I usually do a threshold and compression on the COH to see quality as the analysis proceeds. Or have an array A2(x,y,COH) and as you are using a single sinusoid as the chirp the COH function can be reduced to a single variable.

Also I would rebuild the signal from all my analysed chirps using the position, phase, length and period plus the amplitude and subtract it from the moon to see what is left. Hopefully it will be short impulses which add energy across all frequency bins.

If you add 3 sinusoids into moon then there is going to be 3 peaks. You need to find them all and not just the maximum.

So get the slider code working well, store the arrays A A1 and A2 and then study the results. Look at a T&C plot and get happy with the quality of the analysis, rebuild the signal and compare them. Now you done.

 

[S-LifeLover]

A(x,y,z) where x is the position, y is the length and z is the complex representation of B1 for those parameters.

I do not understand how anyone can build a graph, where one of the axes - a complex number.
Can you give an example of how to build a graph (y=x), where the axis x - complex.
Here is an example where the axis x - not a complex number

x = [1, 2, 3, 4, 5, 6, 7, 8, 9]
y = x;
plot(x, y);

But as I was trying to plot on complex numbers:

x = [1 + 1j, 2 + 1j, 3 + 1j, 4 + 1j, 5 + 1j, 6 + 1j, 7 + 1j, 8 + 1j, 9 + 1j]
y = x;
plot(x, y);

Matlab showed warning:

Warning: Imaginary parts of complex X and/or Y arguments ignored

A (x, y, z) - It seems that this is a 4-dimensional space. With all respect, I do 3-dimensional graph is difficult to build on the 2-dimensional screen.

x = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
y = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10];

A = [];
for xi =1:length(x)
    for yi =1:length(y)
        A(xi, yi) = x(xi) + y(yi);
    end;
end;

mesh(A);

 

[included]

I now have a copy of Matlab reinstalled on my machines so it may be easier in the future to understand what you show in your code. But I do not use Matlab as it cant do stuff without lots of intermediate code. One would like a program to plot the Re or Im or Mag or Phase because it is designated as complex. Whatever. So if Matlab does not do it inherantly then you need to to supply the intermediate code. Its your project after all.

You have to get the slider code right and understand what it does or you cant move forward.

So we have an array A(x,y,z) and x is the position and y is the length.

Now run a loop for differant periods and phases for each of these x and y values and store the answer when it reaches a maximum. That means you change the period and the phase till you get a maxima. Now store the period and phase (in array A2 for later use) if that suites you better.

Of course if you do store the data as complex form it reduces the number of arrays but you have to do more calculations before you can display it - like calculate the magnitude of the complex vector and put it in another array before you display it - as a basic mag function.

Do you understand T&C algorithms and what they show you??? As that may help you in the long run.

This is not a simple task so I would expect the final code to be between 3000 and 6000 lines of a language suited to signal processing. So we have a long way to go yet.