You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Given an audio file, how can I construct a histogram of all bipolar pulses?
2 views (last 30 days)
Show older comments
Hi! I'm looking to take an audio file (.flac) and append each instance of a bipolar pulse to a histogram vector. Is there a version of the "findpeaks" function that only looks for bipolar signals (ie: looks for peaks that are adjacent to a negative peak)? If not, how would I go about constructing a signal template and processing through my file to search for signals of a similar shape? I'm not sure what key words I should be Googling (maybe wavelet or templating?) so any help is appreciated. Thanks!
Accepted Answer
Star Strider
on 14 Nov 2023
The findpeaks function can easily detect negative peaks.
Simply negate the original signal to get them —
t = linspace(0, 50, 125);
s = tan(2*pi*t);
[pks,plocs] = findpeaks(s, 'MinPeakProminence',10); % Peaks
[vys,vlocs] = findpeaks(-s, 'MinPeakProminence',10); % Valleys
Pulse = plocs(plocs-vlocs <= 10)+[-5; 10]; % Pulses
figure
plot(t, s)
hold on
plot(t(plocs), pks, '^r')
plot(t(vlocs), -vys, 'vg')
plot(t(Pulse), [10; 10], '-k', 'LineWidth',2)
hold off
grid
axis('padded')
I do not understand the histogram reference.
Make appropriate changes to get the result you want.
.
11 Comments
AZ0
on 15 Nov 2023
Moved: Star Strider
on 15 Nov 2023
Hi Star Strider,
Thank you so much for your response! I really appreciate it.
To rephrase what I'm trying to do, I guess I'm looking to get a
[counts, z] = hist(pulse, binsize)
so that I am able to plot the number of bipolar pulses as a function of their pulse height.
Do you think you could explain what you are doing with this line:
Pulse = plocs(plocs-vlocs <= 10)+[-5; 10]; % Pulses
Again, thank you for your help!
Star Strider
on 15 Nov 2023
My pleasure!
First, this line:
Pulse = plocs(plocs-vlocs <= 10)+[-5; 10]; % Pulses
defines each pulse (according to my criteria) and is then used to plot a horizonotal balck line identifying them. It identifies the beginning of a pulse if the difference between the peak and valley are less than 10 index units. Once you have identified the pulses, you could then use the maximum in that signal range (or the peak-to-peak value) to create a vector for the histogram counts function to use.
That might go something like this —
t = linspace(0, 1500, 1250);
s = tan(2*pi*t);
[pks,plocs] = findpeaks(s, 'MinPeakProminence',50); % Peaks
[vys,vlocs] = findpeaks(-s, 'MinPeakProminence',50); % Valleys
Pulse = plocs(plocs-vlocs <= 10)+[-5; 10] % Pulses
Pulse = 2×12
47 149 256 358 460 567 669 776 878 980 1087 1189
62 164 271 373 475 582 684 791 893 995 1102 1204
for k = 1:size(Pulse,2)
Pulsev{k,:} = Pulse(1,k) : Pulse(2,k);
maxPulse(k) = max(s(Pulsev{k}));
p2pPulse(k) = max(s(Pulsev{k})) - min(s(Pulsev{k}));
end
maxPulse
maxPulse = 1×12
159.0255 53.0029 795.1377 72.2807 37.8549 113.5882 46.7657 265.0448 61.1590 34.5616 88.3449 41.8414
p2pPulse
p2pPulse = 1×12
200.8669 141.3478 829.6993 133.4397 302.8997 160.3539 160.3539 302.8997 133.4397 829.6993 141.3478 200.8669
figure
plot(t, s)
hold on
plot(t(plocs), pks, '^r')
plot(t(vlocs), -vys, 'vg')
plot(t(Pulse), [1; 1]*50, '-k', 'LineWidth',2)
hold off
grid
axis('padded')
xlabel('Time')
ylabel('Amplitude (mV)')
figure
histogram(maxPulse, 10)
xlabel('Pulse Maximum')
ylabel('Counts')
title('Pulse Height')
figure
histogram(p2pPulse, 10)
xlabel('Pulse Magnitude')
ylabel('Counts')
title('Pulse Peak-To-Peak')
Thjat seems to work acceptably well with my synthetic data. You will likely have to adjust the 'MinPeakProminence' values (individually for the peaks and valleys) to fit your actual data
.
AZ0
on 24 Nov 2023
Again, Star Strider, thank you so much! Sorry about the late response but I just had a follow up question.
I've been getting the following error:
Arrays have incompatible sizes for this operation.
Error in MatlabForumBipolar (line 19)
Pulse = plocs(plocs-vlocs <= 10);
And I think it's because my signal is not one that's always bipolar. I have more peaks than valleys so I have more plocs than vlocs as a result. This makes sense to me. I'm guessing that if I wanted a vector that works for this histogram, I'd have to iterate through plocs and vlocs for peaks and valleys that are of the same magnitude and at around the same time (not just around the same time as you have done? correct me if I'm wrong). I just wanted to make sure of that before I go down the wrong path. Thanks! And if you celebrate, hope you've had a nice Thanksgiving!
Star Strider
on 24 Nov 2023
It would be helpful to have the signal (at least a representative sample) to work with.
I did my best to simulate the situation you describe, setting the 'MinPeakProminence' in the ‘Valleys’ detection to produce fewer of them than there are ‘Peaks’.
This works with my simulated data —
t = linspace(0, 1500, 1250).';
s = tan(2*pi*t);
[pks,plocs] = findpeaks(s, 'MinPeakProminence',50); % Peaks
[vys,vlocs] = findpeaks(-s, 'MinPeakProminence',125); % Valleys
for k = 1:numel(vlocs)
[dist,ix] = min(abs(plocs-vlocs(k)));
Pulse(k,:) = [NaN NaN];
if dist <= 10
Pulse(k,:) = vlocs(k)+ix+[-5; 10]; % Pulses
end
end
Pulse
Pulse = 5×2
156 171
470 485
681 696
995 1010
1206 1221
for k = 1:size(Pulse,1)
Pulsev{k,:} = Pulse(k,1) : Pulse(k,2);
maxPulse(k) = max(s(Pulsev{k}));
p2pPulse(k) = max(s(Pulsev{k})) - min(s(Pulsev{k}));
end
maxPulse
maxPulse = 1×5
1.6511 1.6234 1.7592 1.7287 1.8783
p2pPulse
p2pPulse = 1×5
89.9960 266.6681 27.3958 33.5237 16.8587
figure
hp{1} = plot(t, s, 'DisplayName','Signal');
hold on
hp{2} = plot(t(plocs), pks, '^r', 'DisplayName','Peaks');
hp{3} = plot(t(vlocs), -vys, 'vg', 'DisplayName','Valleys');
hp{4} = plot(t(Pulse)+[-10 10], [1 1]*200, '-k', 'LineWidth',5, 'DisplayName','Identified Biploar Pulses');
hold off
grid
axis('padded')
xlabel('Time')
ylabel('Amplitude (mV)')
legend([hp{1},hp{2},hp{3},hp{4}(1)],'Location','best')
figure
histogram(maxPulse, 10)
xlabel('Pulse Maximum')
ylabel('Counts')
title('Pulse Height')
figure
histogram(p2pPulse, 10)
xlabel('Pulse Magnitude')
ylabel('Counts')
title('Pulse Peak-To-Peak')
This iterates through the ‘Valleys’ (since there are fewer of them), and then uses the minimum peak distance from each to define a bipolar pulse. I included a distance test for ‘Pulse’ and set the entries that do not meet that test to NaN. Set the distance interval (here 10) to whatever you want.
Also, my code here works equally well with row or column vectors for ‘t’ and ‘s’. (I tested it to be sure.)
Thank you! I do, and I did. I hope you did, too!
.
AZ0
on 26 Nov 2023
Hi! Just receiving your message and taking some time to digest it.
For the following line, are you creating a vector called Pulse and appending an entry with index k: [k, the smaller timestamp of the peak and valley]? If so, why are you adding the index to the valley size and what is the [-5;10] for?
Pulse(k,:) = vlocs(k)+ix+[-5; 10];
And could you explain what the for loop iterating through the size of Pulse does?
Sorry about all of the confusion. Also I've attached a very short sample of the signal I'm trying to work with in the following link. Please let me know if you have any trouble downloading it.
Thanks!
Star Strider
on 26 Nov 2023
My pleasure.
I prefer not to download files from external sites because I cannot download them to here and I do not want to download them to my computer and then upload them here. (Ever since the Run utility — the right-facing green arrow in the top toolstrip — was added here, I have used it or MATLAB Online almost exclusively to develop and run code on MATLAB Answers.) Use the ‘paperclip’ icon to the right of the Σ to upload them here. If they have an unrecognised extension (such as .dat), use the zip function to enclose it in a .zip file and then upload the .zip file.
I will wait for the file to appear here.
The ‘Pulse’ loop iterates through ‘vlocs’ because you mentioned that there are fewer valleys than peaks. The ‘Pulse’ assignment begins with the ‘vlocs’ index, adds ‘ix’ to it because the min function returns the index (‘ix’) relative to the beginning of the difference between the ‘plocs’ vector and the indexed ‘vlocs’ value, so adding ‘ix’ produces the absolute index (for the entire vector). The loop counter ‘k’ is the ‘vlocs’ index.
The size-of-Pulse loop uses the ‘[-5; 10]’ offset vector result in ‘Pulse’ to create a range to determine the maximum and peak-to-peak values within that range. (I chose the ‘[-5; 10]’ range because it works with my simulated data. It may need to be changed to work with real data.)
.
AZ0
on 26 Nov 2023
Moved: Star Strider
on 26 Nov 2023
Ok, I see. I think that makes sense to me. Thanks!
I've attached the file below, please let me know if there are any issues with it!
Star Strider
on 26 Nov 2023
My pleasure!
After a bit of tweaking of my original code, this works —
% t = linspace(0, 1500, 1250).';
% s = tan(2*pi*t);
Uz = unzip('matlab sample audio.flac.zip')
Uz = 1×1 cell array
{'matlab sample audio.flac'}
[s, Fs] = audioread(Uz{1})
s = 2128109×2
-0.0085 0.0378
-0.0055 0.0339
-0.0032 0.0305
0.0019 0.0240
0.0073 0.0196
0.0109 0.0227
0.0109 0.0289
0.0156 0.0273
0.0191 0.0294
0.0172 0.0304
Fs = 4410
L = size(s,1);
t = linspace(0, L-1, L).'/Fs;
figure
plot(t, s(:,1), '-', 'DisplayName','Left Channel')
hold on
plot(t, s(:,2), 'DisplayName','Left Channel')
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
title('Signals Plot')
legend('Location','best')
SelectIdx = [6483 6512];
SelectT = t(SelectIdx);
SelectWidth = diff([1556 1606])
SelectWidth = 50
figure
plot(t, s(:,1), '.-', 'DisplayName','Left Channel')
hold on
plot(t, s(:,2), 'DisplayName','Left Channel')
hold off
grid
xlabel('Time (s)')
ylabel('Amplitude (Units)')
xlim(SelectT+[-0.015; 0.015])
xline(SelectT(1), '--c', 'DisplayName','Valley')
xline(SelectT(2), '--g', 'DisplayName','Peak')
legend('Location','best')
title('Arbitrary ‘Zoom’ To Examine Signals Detail')
s = s(:,1); % Choose Left Channel
[pks,plocs] = findpeaks(s, 'MinPeakProminence',0.2); % Peaks
[vys,vlocs] = findpeaks(-s, 'MinPeakProminence',0.2); % Valleys
ofst = [-50 50];
for k = 1:numel(vlocs)
[dist,ix] = min(abs(plocs-vlocs(k)));
Pulse(k,:) = [NaN NaN];
if dist <= 100
idxrng = vlocs(k)+ix+ofst;
if idxrng(1) < 1
idxrng(1) = 1;
end
if idxrng(2) > L
idxrng(2) = L;
end
Pulse(k,:) = idxrng; % Pulses
end
end
NrNaN = nnz(isnan(Pulse(:,1)))
NrNaN = 105
% Pulse
for k = 1:size(Pulse,1)
maxPulse(k,:) = NaN;
p2pPulse(k,:) = NaN;
if ~all(isnan(Pulse(k,:)),2) & ~all(isempty(Pulse(k,:)),2)
Pulsev = Pulse(k,1) : Pulse(k,2);
if numel(Pulsev) < 2 % Check For Empty 'Pulsev'
break
end
maxPulse(k,:) = max(s(Pulsev));
p2pPulse(k,:) = max(s(Pulsev)) - min(s(Pulsev));
end
end
% maxPulse
% p2pPulse
figure
hp{1} = plot(t, s, 'DisplayName','Signal');
hold on
hp{2} = plot(t(plocs), pks, '^r', 'DisplayName','Peaks');
hp{3} = plot(t(vlocs), -vys, 'vg', 'DisplayName','Valleys');
% hp{4} = plot(t(Pulse)+[-10 10], [1 1]*200, '-k', 'LineWidth',5, 'DisplayName','Identified Biploar Pulses');
hold off
grid
axis('padded')
xlabel('Time')
ylabel('Amplitude (mV)')
% legend([hp{1},hp{2},hp{3},hp{4}(1)],'Location','best')
legend([hp{1},hp{2},hp{3}],'Location','best')
figure
histogram(maxPulse, 5E+2, 'EdgeColor','none')
xlabel('Pulse Maximum')
ylabel('Counts')
title('Pulse Height')
figure
histogram(p2pPulse, 5E+2, 'EdgeColor','none')
xlabel('Pulse Magnitude')
ylabel('Counts')
title('Pulse Peak-To-Peak')
Not much needed tweaking, since my original code accounted for most of it. I had to change the ‘ofst’ index range to accommodate the sampling intervals with respect to the peaks, and adjusted the detection limits in the findpeaks calls. (Those use my favourite name-value pair argument —'MinPeakProminence' — however you can choose any criterion and detection level level that works best for you.) Some other changes were necessary with respect to the indexing (it took a few minutes to figure out where that problem — indexing beyond the ranges of the vectors — actually was, fixing it involved an if block to check ‘Pulsev’) and with those tweaks, it works. It should do what you requested.
.
AZ0
on 27 Nov 2023
Got it, thank you so much! That seems really helpful. I'll try to make sense of it.
Star Strider
on 27 Nov 2023
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
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: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- 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)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)