Given an audio file, how can I construct a histogram of all bipolar pulses?

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!

Answer 1

[Star Strider]

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.

.

[AZ0]

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]

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]

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]

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]

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.

https://drive.google.com/file/d/1-xK7m7oLv2wC7KbdbIR4fFXZauu1BJs2/view?usp=sharing

Thanks!

[Star Strider]

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]

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]

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]

Got it, thank you so much! That seems really helpful. I'll try to make sense of it.

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

[Star Strider]

Thank you!