MATLAB Answers

- You will see updates in your activity feed.
- You may receive emails, depending on your notification preferences.

12 views (last 30 days)

Sign in to answer this question.

Star Strider
on 10 Feb 2020

steamrice
on 10 Feb 2020

Thanks for the response! That may work! I am just looking at an example

findchangepts(t,'MaxNumChanges',8)

But should I segment a way such that it seperate, says the 8 events, from the inital point to the end point of each event? (would it be possible?)

Star Strider
on 10 Feb 2020

My pleasure!

Without your signal to experiment with, it is not possible for me to determine what the best option would be. Consider using other name-value pairs as well, such as 'Statistic' and others so that findchangepts will do what you want more easily. There are similar functions that you may consider experimenting with, such as cusum and ischange (R2017b and later), both with different calling conventions and syntax.

Star Strider
on 10 Feb 2020

Using:

D = load('pilotexample.mat');

data = D.data;

EMG = data(:,1);

Ts = 0.5;

t = linspace(0, size(data,1), size(data,1)).'*Ts;

cpv = findchangepts(EMG, 'MaxNumChanges',20,'Statistic','rms', 'MinDistance',fix(numel(EMG)/500))

the findchangepts output for the first column of data are:

cpv =

6408

8492

13240

15902

20901

22658

27312

28504

41208

42055

47777

49485

54219

56000

61036

61889

66853

67899

Elapsed time is 555.593266 seconds.

corresponding to:

I will re-run this with 22 change points (for the other 2 spikes).

steamrice
on 10 Feb 2020

steamrice
on 10 Feb 2020

Star Strider
on 10 Feb 2020

As always, my pleasure!

My apologies for the delay — I had an errand to run.

The complete code (now that I have it working):

D = load('pilotexample.mat');

data = D.data;

EMG = data(:,1);

Ts = 0.5;

t = linspace(0, size(data,1), size(data,1)).'*Ts;

cpv = findchangepts(EMG, 'MaxNumChanges',24,'Statistic','rms', 'MinDistance',fix(numel(EMG)/500))

tic

EMG = data(:,1);

Ts = 0.5;

t = linspace(0, size(data,1), size(data,1)).'*Ts;

cpv = findchangepts(EMG, 'MaxNumChanges',24,'Statistic','rms', 'MinDistance',fix(numel(EMG)/500))

toc

figure

plot(t, data(:,1))

hold on

yl = ylim;

plot(t([cpv cpv].'), (ones(size(cpv))*ylim).', 'r')

hold off

grid

for k = 1:fix(numel(cpv)/2)

segment{k} = [t(cpv(k):cpv(k+1)), EMG(cpv(k):cpv(k+1))]; % Save Segments In Cell Array

end

figure

subplot(2,1,1)

plot(segment{1}(:,1), segment{1}(:,2))

grid

title('Segment 1')

subplot(2,1,2)

plot(segment{10}(:,1), segment{10}(:,2))

grid

title('Segment 10')

with the segmented data plot:

and:

This gives you a way of saving the segments and then addressing them as I do in figure(2). The code creating figure(1) plots the change points as vertical red lines.

Yes, you would have to use findpeaks on each ‘segment’ (as I call them here) to get the maximia. You might be able to do that in the loop that creates the individual ‘segment’ cell array. (I suggest saving the ‘segment’ array to a .mat file if you want to use it again. It would save the time of re-segmenting it.)

It takes half as long with 24 change points as with 20, likely because findchangepts does not have to ‘decide’ which 20 to report in its output. Using 24 also gets all of them. (Unfortunately, it may not be possible to completely automate a way to correctly determining the number of change points.)

steamrice
on 11 Feb 2020

Star Strider
on 11 Feb 2020

As always, my pleasure!

They are representing the absolute times. (I opted for that in order to preserve them in the event that was necessary.) If you want to represent the relative times (so that each starts at 0), subtract the initial time from the rest of the vector.

For example:

plot(segment{10}(:,1)-segment({10}(1,1), segment{10}(:,2))

will plot ‘segment{10}’ beginnining at 0. That applies to all the others. Using relative times might be easier if you are using findpeaks. It would also standardise them if you are interested in exploring morphological similarities.

steamrice
on 11 Feb 2020

That helps! I was just testing and learning how to find peaks for the next steps. But I am a bit confused. I tried : (say I wanted to find peak of the first segment)

findpeaks(segment{1) (:,2),'MinPeakHeight',0.1 )

The plots showed the segment 1 and give peak of the amplidue that is higher than 0.1. So for segment 1, I should use segment{1) (:,2) right instead of segment{1) (:,1)? I am a bit confused here.

Also, what did you by "You might be able to do that in the loop that creates the individual ‘segment’ cell array." for the findpeaks?

Currently I have to kinda plot the graph to see its ampliude to set the MinPeakHeight. Is it the standard way to do it? I am wondering would there be a more efficient way to find that amplidue if the signals are differnet

Star Strider
on 11 Feb 2020

I was thinking of something like this:

for k = 1:fix(numel(cpv)/2)

segment{k} = [t(cpv(k):cpv(k+1)), EMG(cpv(k):cpv(k+1))]; % Save Segments In Cell Array

[pks{k},locs{k}] = findpeaks(EMG(cpv(k):cpv(k+1)), 'MinPeakHeight',sqrt(mean(EMG(cpv(k):cpv(k+1)).^2)));

end

figure

subplot(2,1,1)

plot(segment{1}(:,1), segment{1}(:,2))

hold on

plot(segment{1}(locs{1}), pks{1}, '^r')

hold off

grid

title('Segment 1')

subplot(2,1,2)

plot(segment{10}(:,1), segment{10}(:,2))

hold on

plot(segment{10}(locs{10}), pks{10}, '^r')

hold off

grid

title('Segment 10')

doing each findpeaks call as the ‘segment’ arrays are created, then producing this plot:

One of the peaks in ‘segment{1}’ definitely appears to be higher than 0.1. Is that a problem? (Note that if there is a strict limit on the amplitudes that you expect, the apparent ‘overshoot’ could be due to a sampling artifact. In my own studies, a value that occurs at the same time as a sampling time in a rapidly-changing signal can be distorted by the sampling process. I have no idea why that occurs, I only know that it does.)

You can extend the subplot figure to include all of the ‘segment’ vectors. Since there are 12 in this record, having 6 rows of subplots and 2 columns might work.

I used the RMS value of each EMG ‘segment’ vector for 'MinPeakHeight', since the amplitudes vary significantly in the EMG between the ‘segment’ vectors, and that appeared to be a relatively reasonable and adaptive approach. Use multiples of the RMS value to get different results, or your own value for 'MinPeakHeight' and other name-value pairs that findpeaks offers. The amplitudes appear to having been returned correctly here.

steamrice
on 11 Feb 2020

Thanks for the fast response! I will playing around and give you some updates!

Star Strider
on 11 Feb 2020

As always, my pleasure!

I wish it could be faster, however findchangepts — for all its wonderful abilities — still takes 4½ minutes for this data set.

steamrice
on 12 Feb 2020

Thanks for the response! I think this way would better help me in terms of getting the amplitude because potentially the signal could vary signifcantly. Having manually set the ampltude (like 0.1 for segement would not be efficient!)

I was looking back to the 'segement', for this line:

cpv = findchangepts(EMG, 'MaxNumChanges',20,'Statistic','rms', 'MinDistance',fix(numel(EMG)/500))

Q1: I am not too sure how did you come up with fix(numel(EMG)/500). Why did you divide by 500?

Q2: Could I also think of those segment as 'onset' of the EMG signals (as the stimulus occurs)?

Star Strider
on 12 Feb 2020

‘I am not too sure how did you come up with fix(numel(EMG)/500). Why did you divide by 500?’

I adapted it to the length of the record. Use any value you want.

‘Could I also think of those segment as 'onset' of the EMG signals (as the stimulus occurs)?’

Yes. The first change point in each segment would work for that purpose. (I initially thought the trigger channel (third column) might be useful. However the third column of your .mat file only appeared to have two triggering events, so I did not use it. Also, they did not correspond to the other data.)

In terms of speed, I would use 24 change points for this record. That takes half the time that 20 change points does.

I just now thought of a way to automate determining the number of segments (half the number of change points):

EMG = data(:,1);

Ts = 0.5;

t = linspace(0, size(data,1), size(data,1)).'*Ts;

env = envelope(EMG, 2500, 'rms');

[epks,elocs] = findpeaks(env, 'MinPeakProminence', 0.0005);

plot(t, EMG)

hold on

plot(t, env)

plot(t(elocs), epks, '^g')

hold off

grid

<REST OF THE CODE>

If you have the Signal Processing Toolbox (R2015b and later), you have the envelope function.

For this record, it returns 12 peaks, so it does count the correct number of EMG events. Experiment with that to see if it works correctly with your other records. It might be necessary to change the 'MinPeakProminence' value, however the 0.0005 value works with this record.

steamrice
on 12 Feb 2020

Star Strider
on 12 Feb 2020

As always, my pleasure!

I wish I had thought of it sooner! It is still necessary to use findchangepts, however with the correct number of change points for each record (2*numel(epks)) that should be easier.

steamrice
on 12 Feb 2020

You have been so helpful!! I am just looking into more detail regarding the segment related to the number of the events.

(I initially thought the trigger channel (third column) might be useful. However the third column of your .mat file only appeared to have two triggering events, so I did not use it. Also, they did not correspond to the other data.)

The trigger signal that I sent was supposed to correlate when a stimuls happens. It's not so useful at the moment as I pogrammed it only 'logic high' when the stimuls happens 'during that stimlus period' using Matlab. So if you see the first picture I posted you would see there is first 4 events and pause and the last 4 events (Total of 8 events). The trigger signal will be high for the first 4 and low during the pause and high again for the last 4. This is because I put the trigger siganl within the events in Matlab. I am hoping to change it in such a way that the trigger signal only high when there is an event from the EMG. (Like a series of logic high and logic low for each event instead of each 'period'). I do not know whether it will be possbile....(This may affect the way I analyze the data).

For

env = envelope(EMG, 2500, 'rms');

The sliding window of length wl samples (2500) in this case. How should we determine this number? Is it dependent on the frequency of the signal? (I think it was 2000Hz)

Thanks again!

Star Strider
on 12 Feb 2020

As always, my pleasure!

The window length for envelope is what I determined experimentally to provide a reasonable result. Like most everything in signal processing with physiological signals, the parameter choices are chosen by experimenting and then using the parameter that appears to give the best result. There is nothing inherently ‘magic’ about 2500, so experiment with other values. You may find one that works better.

steamrice
on 12 Feb 2020

Thanks for your time! Do you have suggestion for the trigger siganls from your pervious experience?

Star Strider
on 12 Feb 2020

In the one multiple-animal study that I did, the stimuli were immediately followed by a muscle contraction, so the stimuli (recorded on a different channel) could be used as the start of the subsequent EMG record, and the end of the previous record. That made segmenting the EMG signals much easier. That does not appear to be the same for your data. Also, I thought ‘trigger’ = ‘stimulus’. It does not with your data.

No worries about that — you designed your study diifferently, and apparently with a different objective, than I did mine.

steamrice
on 12 Feb 2020

Star Strider
on 12 Feb 2020

steamrice
on 13 Feb 2020

WOW. That's impressive in the 1990s!

That might be a liitle be unrelated to what I orginally asked, but I am getting struck and hope to get some ideas from you. For the EMG measurement, I have an external device connected to the matlab that can send the logic high and logic low as a marker corresponding to the stimulus. But I think it does not work too well. What I was thinking was hoping to get the trigger signal as event marker to tell when the stimulus happens and when the event was being carried out. (event happens as logic high and stimulus being carried out as logic low). This is because the elasped time running for the Matlab code is not matching the stimulus. When you were woring with EMG measruemnt, did you have a way/were you interested in finding/ getting the correct stimulus repsonese's time ?

Star Strider
on 13 Feb 2020

steamrice
on 13 Feb 2020

Star Strider
on 13 Feb 2020

I believe so, however the ‘trigger’ does not seem to be an actual stimulus. That may be the difference.

I was simply describing the setup I used.

steamrice
on 13 Feb 2020

That's what I was thinking and planning to do. Is that something closer to what you meant?

Star Strider
on 13 Feb 2020

steamrice
on 14 Feb 2020

Star Strider
on 14 Feb 2020

steamrice
on 14 Feb 2020

I see! I am playing things around again and will ask you for more questions! Thank you!

steamrice
on 18 Feb 2020

Good day!

I was playing around my setup and realized that the problem I have now is more related to the timing setup! I am worried about is the Matlab script that I have for running the experiment. I have a "for loop" that has:

getting the audio file>saving the audio file (by calling a function)>playing the audio (by calling another function). I'd like to have a 'trigger signal' such that the stimulus plays (1s) and record the onset of the person (1s).

But I just checked the timing of the 'for loop' I have by using tic and toc. It says 'Elapsed time' is about 3.3/3.4s for the loop to run. I am thinking that the unmatching timing comes from the MATLAB script... I am unsure of how to fix this at this point since each command seems necessary (I want to keep what I have: opening the audio, saving the audio, playing the audio with defined frequency...)

I have attached a picture trying to show you what my problem is.

Thank you again for your time!

Star Strider
on 18 Feb 2020

steamrice
on 4 Mar 2020

steamrice
on 4 Mar 2020

Star Strider
on 4 Mar 2020

steamrice
on 4 Mar 2020

Thanks for your repsosne as always!

Baesd on the definiiton: "In classifying speech sounds, the lowest strong frequency band is termed the first formant. The next highest is termed the second formant, and so on."

I have thought about using spectrogram as well to have a plot in terms of frequency vs time. But I am not too sure how to detect the frequency that matchs to what I am looking for.

what I am thinking and having is to load the audio, then plot it in terms of dB vs freuency. Then I am just wondering how to select the appropriate "formant". This is what I am trying to acheive:

Creadits to:

I have read a post and thats what I have:

[audio_file, Fs] = audioread("testing.wav");

L = 10000;

Ts = 1/Fs;

Fn = Fs/2;

FT_af = fft(audio_file)/L;

Fv = linspace(0, 1, fix(L/2)+1)*Fn;

Iv = 1:numel(Fv);

[PksL,LocsL] = findpeaks(abs(FT_af(Iv,1))*2, 'MinPeakHeight',0.006, 'MinPeakDistance',20);

[PksR,LocsR] = findpeaks(abs(FT_af(Iv,2))*2, 'MinPeakHeight',0.006, 'MinPeakDistance',20);

figure(1);

semilogx(Fv, 20*log10(abs(FT_af(Iv,2))*2)), xlabel('Frequency, Hz'),ylabel('magnitude');

hold on

plot(Fv(LocsR), PksR, '^r', 'MarkerFaceColor','r')

hold off

xlabel('Frequency (Hz)')

ylabel('Amplitude')

However, the scale is not appropriate yet and I haven't got the the formant properly..

Star Strider
on 4 Mar 2020

‘But I am not too sure how to detect the frequency that matchs to what I am looking for.’

Nor am I. This is quite remote from my areas of expertise.

I would solve it by creating different filters (the firls function could be appropriate here), and then take the time-integrated output from all the filters that filtered the same signal and look for the greatest result for a specific time range of a signal. That result would match the filter for that formant.

Wavelets could be another way to solve this, however I am not sufficiently familiar with wavelets to craft a specific approach.

steamrice
on 4 Mar 2020

Thanks for your resposne!!

I am wondering that is what is the y-axis represents for when we plot the audio in fft? We do have to use

semilogx(Fv, 20*log10(abs(FT_af(Iv,2))*2))

to convert it to dB right?

I would solve it by creating different filters (the firls function could be appropriate here), and then take the time-integrated output from all the filters that filtered the same signal and look for the greatest result for a specific time range of a signal. That result would match the filter for that formant.

Could you direct me on this? Creating different filiters (like low pass and high pass ) and select the frequency range that we are interested in using firls? I am unsure the what do you mean by take the time-integrated output ...

Since I just want to compute the first two formants, will findpeak do the job by looking for the first two lowest strong band in the dB vs Hz plot?

Sorry for asking many questions. Your suggestions have really pushed me to exlore different options and helped me to formule possible solutions!

Star Strider
on 4 Mar 2020

As always, my pleasure!

Doing this transformation:

20*log10(abs(FT_af(Iv,2))*2)

converts ‘FT_af’ to dB.

See the documentation for the firls function to understand how to create the filters. After you design them (using firls or designfilt), use freqz to be sure the filter works as you intend it to work.

The time integration simply uses trapz (or cumtrapz depending on how you formulate the problem) to integrate the outputs of the various filters. The filter with the greatest output is likely the one you want. As I mentioned, this is very far from my areas of expertise, so there may be much more efficient ways of doing this.

If you want to use findpeaks, it will be necessary for you to take the Fourier transforms of specific regions of the time-domain vector that represent the sounds you want to examine. Taking the Fourier transform of the entire record will lllikely give you no useful information.

steamrice
on 4 Mar 2020

Star Strider
on 4 Mar 2020

Plot the peaks in dB:

plot(Fv(LocsR), 20*log10(PksR), '^r', 'MarkerFaceColor','r')

They should now be on the curve.

Also, for your purposes, it might be best to use a linear frequency axis rather than a logarithmic one. It may make the information easier to see.

steamrice
on 4 Mar 2020

Also, for your purposes, it might be best to use a linear frequency axis rather than a logarithmic one. It may make the information easier to see.

Do you mean fft(PksR) in my case?

Star Strider
on 4 Mar 2020

No. The ‘PksR’ vector are the peaks of the absolute value of the fft, so just plot them using the decibel transformation. Nothing else is necessary.

The linear frequency axis is simply a suggestion in order to make the spectrum easier to visualise.

steamrice
on 4 Mar 2020

steamrice
on 4 Mar 2020

I think I am still trying to find out a way to better do this. If I just want to find the first two higest peaks from the fft plo in a certain freqency range.(does that make much sense to you?) Is it possbile to just use findpeaks() to do so? You taguht me how find peaks worked and I tried to use in this case

I am thinking:

[audio_file, Fs] = audioread("testing.wav");

L = 10000;

Ts = 1/Fs;

Fn = Fs/2;

FT_af = fft(audio_file)/L;

Fv = linspace(0, 1, fix(L/2)+1)*Fn;

Iv = 1:numel(Fv);

[PksR,LocsR] = findpeaks(abs(FT_af (:,1)), 'MinPeakHeight',0.006, 'MinPeakDistance',20);

plot(Fv(LocsR), abs(FT_af (:,1)), '^r', 'MarkerFaceColor','r');

But I recived the following error:

Index exceeds the number of array elements (5001).

Error in

plot(Fv(LocsR), abs(FT_af(:,1)), '^r', 'MarkerFaceColor','r');

Star Strider
on 4 Mar 2020

Using findpeaks makes sense. I am not certain what peaks you are looking for.

Experiment with 'MinPeakProminence' to detect them. It may be more reliable than 'MinPeakHeight' and 'MinPeakDistance'.

steamrice
on 4 Mar 2020

Another thing if you don't mind, I was trying the Ipc method:

load testing

segmentlen = 100;

noverlap = 90;

NFFT = 128;

spectrogram(FT_af (:,1),segmentlen,noverlap,NFFT,Fs,'yaxis')

title('Signal Spectrogram')

dt = 1/Fs;

I0 = round(0.1/dt);

Iend = round(0.25/dt);

x = (I0:Iend);

x1 = x.*hamming(length(x));

preemph = [1 0.63];

x1 = filter(1,preemph,x1);

A = lpc(x1,6);

rts = roots(A);

rts = rts(imag(rts)>=0);

angz = atan2(imag(rts),real(rts));

[frqs,indices] = sort(angz.*(Fs/(2*pi)));

bw = -1/2*(Fs/(2*pi))*log(abs(rts(indices)));

nn = 1;

for kk = 1:length(frqs)

if (frqs(kk) > 90 && bw(kk) <400)

formants(nn) = frqs(kk);

nn = nn+1;

end

end

formants

I am having a trobule with A = lpc(x1,6);

Error using roots (line 23)

Input must be a vector.

Error in testingmethod(line 77)

rts = roots(A);

Star Strider
on 4 Mar 2020

Sign in to comment.

Sign in to answer this question.

Opportunities for recent engineering grads.

Apply TodayUnable to complete the action because of changes made to the page. Reload the page to see its updated state.

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

Select web siteYou can also select a web site from the following list:

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

- América Latina (Español)
- Canada (English)
- United States (English)

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

## 0 Comments

Sign in to comment.