Extracting multiple envelopes from a signal

I have a signal which when observed carefully shows a three individual sinusoidal signals (figures attached). I want these three individual signals to be extracted in terms of amplitude, frequency and phase. Please help!

[John BG]

would you please be so kind to supply the signal?

in a .mat file

thanks in advance, awaiting answer

[First Last]

Sure! Here is the .mat file. Thanks!

Answer 1

[Cedric Wannaz]

Quick cheap trick gets 2/3. Was going to bed, no time for more but that could be a starting point.

 function main

    load( 'MultipleSignalEnvelope', 'ITDc0' ) ;

    plot( ITDc0, 'b.' ) ;
    grid on ;
    hold on ;

    t = 1 : numel( ITDc0 ) ;
    A = 1.2 ;

    P = fminsearch( @objective, zeros(3,1) ) ;

    plot( A*sin((t+P)/360*2*pi).' ) ;

    function mindif = objective( P )
        mindif = sum( min( abs( ITDc0-A*sin((t+P)./360*2*pi )))) ;
    end

 end

%

[First Last]

This is nice, but you assumed constant amplitude of 1.2 for all three waveforms. This may not be always true. These waveforms might have different amplitudes. Also, the number of waveforms may vary.

[Cedric Wannaz]

Then you should try to re-post your question using your real name, reference this thread, and give much more explanation initially to give people more incentive to spend time on your question.

The general case is not a trivial task and sadly we get a lot of students who come anonymously with their semester project and ask "solve the general case for me" without having spent much time trying/learning, and thinking about the complexity of the general case (e.g. I have to read license plates from all countries in all weather conditions and by night, please provide code).

I'm not saying that you'll get an answer, but if you want to increase your chances that someone with experience in signal processing will spend time on your thread, don't make your question look like the example above.

[First Last]

This is not a semester project. But I guess, you are correct. I'll go ahead and repost this with details. Thanks for the attempt.

[Cedric Wannaz]

No problem. The issue was not really you here but more the background noise on the forum. Many students create anonymous accounts quickly, ask for the world and then disappear. This doesn't encourage us to spend a lot of time on complex issues (hence ImageAnalyst answer I guess) unless the person and the question look serious.

[First Last]

Thanks Cedric, I got it. Just a quick question. The code above gives an error as t and P cannot be added due to different dimensionality.

[Cedric Wannaz]

You are probably using a version of MATLAB older than 2016b. From 2016b on, MATLAB supports automatic expansion. We can check this quickly by evaluating:

 >> [1;2;3] .* [10,20,30,40]
 ans =
    10    20    30    40
    20    40    60    80
    30    60    90   120

On your computer it should generate an error(?) What happens here is that MATLAB expands automatically the operands "along the dimension of each other" in a fashion that allows the operation. I'm not wording that too well, but you can understand looking at the output, which is what you would get with:

 1 1 1 1               10 20 30 40
 2 2 2 2      .*       10 20 30 40 
 3 3 3 3               10 20 30 40

On former versions of MATLAB, we had/have to use BSXFUN for this ( an explicit call to a function that does an implicit expansion ;) ). I cannot test, but if you replace t+p by bsxfun( @plus, t, P ), that should solve a first expansion issue, and then same thing for ITDc0-A.... So overall you should replace:

 ITDc0-A*sin((t+P)./360*2*pi)

by

 bsxfun( @minus, ITDc0, A*sin( bsxfun( @plus, t, P)./360*2*pi )

Given my first explanation and what you understand of the approach, you probably understand that P is a vertical vector of three "phases" (term that I am loosing loosely here because there are multiplicative factors), t+P (now built using BSXFUN) is an array of three horizontal time vectors stacked vectically that will "de-phase" the sine, ITDc0-A.. a vector of differences between your signal and the three sines, and the FMIN call is aiming at finding the three phases that minimizes the min absolute difference between the signal and the three sines.

Now this is not a good approach; it has nothing to do with digital signal processing, it is just minimizing something that is not even that smart because it fails 1/3. But it was a fun attempt for the 5 minutes that I had.

[First Last]

Thanks! It works now :-)

[First Last]

Do you think we can use the same logic for Amplitude too?

[Cedric Wannaz]

You will have to try. Define P as a vector with 6 components instead of 3, the first three being amplitudes and the next three being "phases". In the objective function, replace A with P(1:3) and P with P(4:6). Also, update the initial conditions. Instead of zeros(3,1), try to start e.g.. with three unit amplitudes and three null phases: [ones(3,1);zeros(3,1)].

Note that you could use FMINCON and define constraints, e.g. upper and lower bounds for the amplitude and phases.

Also note that you'd better try to see if you could get help from someone with experience in digital signal processing, because this approach is just a quick trick. If you reword your question in a new thread and nobody answers, you may be able to hire someone for not that much on sites that market freelancers.

Answer 2

[Image Analyst]

Look at the fft.

[First Last]

fft does not work. I tried. Can you please try and let me know if it does?

[Image Analyst]

Let's see what you tried.

[First Last]

I tried this. Let me know if I am doing something wrong. Thanks!

X=fft(signal);
N = length(signal);
[i,j]=max(abs(fftshift(fft(X))));
amp=i;
freq=wn(j);
amplitude=amp*2/N;

%%extracting the phase
Xp = 1/N * fftshift(fft(signal));
P = atan2(imag(Xp), real(Xp)) * 180/pi;  % Too noisy to see the peaks
Xp(abs(Xp) < 0.0001) = 0;        % Reduce noise
P2 = atan2(imag(Xp), real(Xp)) * 180/pi;
ph=max(P2)*pi/180;

t = (0:1:length(signal)-1);

figure
plot(signal,'r')
hold on
plot(t,amplitude*sin(freq*t+ph),'b')

[Image Analyst]

How did the signal get so messed up in the first place? Are you willing to have some manual interaction to get the job done, or is it "either automatic or forget it"?

[First Last]

Any clue?

[Image Analyst]

Are you going to answer my questions first?

[Image Analyst]

Can you then simply plot it and call ginput(2) to ask the user to visually identify and click on the two peaks that mark the end of a particular cycle/period/wavelength? Once you know the wavelength for all 3 waveforms, then I think you're done with this part and you continue on with the rest of your program.

[First Last]

There are no user inputs allowed. It should be all autonomous.

[Image Analyst]

Well that's a change from what you said before. Too bad, that would have been so easy. I don't know how to solve it without spending way more time on it than I allow myself to spend on questions here, so . . . good luck.

[First Last]

I am sorry, I miss understood your earlier comments. It should be autonomous.