Code for signal analysis

I trying to analyse some electric signal data but I'm still a Matlab beginner and need help. I use the code bellow just to plot time and signal, but i want to stract the repetition rate of the signal (Hz/number of pulses per second). I think this is possible if i delimite a threshold for the start. Other problem is: i have two signals together in the same file (one big, other little, a example in the figure). Its possible to calculate the repetition rate for each signal? Follow a example of my data.

Thanks for reading

file = 'filename';
signal=wavread(file); % signal
fs = 48000; % sample rate 
t=[1/fs:1/fs:length(signal)/fs];
plot(t,signal);

[John BG]

Thiago

You are not supplying enough information.

the signal.mat is not a signal, but 2, without time reference, without the time reference you keep the readers doing guesswork, let me explain:

s=load('signal.mat');
y1=s.signal(:,1);y2=s.signal(:,2);Y=[y1 y2];

now listen to them

Fs=40000;signalObj_Y1=audioplayer(y1,Fs);play(signalObj_Y1)
Fs=55000;signalObj_Y1=audioplayer(y1,Fs);play(signalObj_Y1)
Fs=155000;signalObj_Y1=audioplayer(y1,Fs);play(signalObj_Y1)

As you can hear, attempting to find the Pulse Repetition Frequencies (PRF1 and PRF2) without timebase, is pointless.

Then one may consider that from the graph above:

 10 pulses (the larger amplitude signal), 8 intervals,  of 50ms, this makes PRF1 ~ 25Hz

 14 pulses of the smaller amplitude signal, on same really short 8 intervals of 50ms each, this is, PRF2 ~ 35Hz

So, where are the 50Hz (Euroland) or 60Hz (US) ?

If you are analysing conducted signals along an electric line you should take a frequency scan well above the audio limits, because just with ADSL and car engine spark plugs, and the whole lot of mobile phone towers around, if not any wired Ethernet plugged to a nearby mains, you really have a lot of noise and interference that is not present in the above graph,

So could you please be so kind to supply the time base of the audio signal?

Question: why do you supply 2 signals in the signal.mat?

Observation:

Y_fft1=fft(y1);plot(abs(Y_fft1));

the short graph with 50ms intervals looks like the signal with slower PRF, this is PRF1~25Hz is has stronger amplitude, yet out of the fft of the complete signal supplied in signal.mat it appears to be the other way around, doesn't it?

So the signal with 2 clear tones, or buch of tones, seems to have both wide variations of amplitude and frequency along the supplied sample, would it make sense to also ask for f1_min f1_max f2_min and f2_max?

Awaiting answer

John BG <mailto:[e-mail] [e-mail]>

[Thiago Petersen]

Hi John, Thanks for your help. So, my main objective with this data is just know the mean rate of the respectives spikes (pulses/seconds, S1 = tall spikes; S2 = little spikes). But i'm also trying to find other variables in the data, like the max/minimum frequency (just time based variables). My time refference its the sample rate of the recorder (48000 Hz), aprox 2 seconds in the attached file. Yeah, first i recorded 2 signals in the file, but follows attached a mono recorded file. Thanks!

Answer 1

[Image Analyst]

I'd take the absolute value of the signal (to flip the bottom part up), then I'd threshold at 0.3 and 0.06 to get the two parts of the signal - the high spikes and the shorter spikes. Then use bwlabel to give an ID number to each spike and use regionprops to compute the centroid of each thresholded spike. Once you know that, you can get the average time between spikes, or whatever other information you might want. Untested code (because you didn't include your data) is below:

bigSpikes = abs(signal) > 0.3;
[labeledSpikes, numberOfSpikes] = bwlabel(bigSpikes);
props = regionprops(labeledSpikes, 'Centroid');
allCentroids = [props.Centroid];
xCentroids = allCentroids(1:2:end)
meanSpacingTall = diff(xCentroids)

Do the same for shortSpikes

shortSpikes = abs(signal) > 0.06;
[labeledSpikes, numberOfSpikes] = bwlabel(shortSpikes);
props = regionprops(labeledSpikes, 'Centroid');
allCentroids = [props.Centroid];
xCentroids = allCentroids(1:2:end)
meanSpacingShort = diff(xCentroids)

[Thiago Petersen]

Hi Image Analyst! Thanks for your help! You knows a lot. I'm trying to make it work here. Your code partially worked, but i have some troubles to identify some big/little spikes. To you have a better idea of my data, i attached a part of the signal in this comment. Please can take a look if this code is able to identify the spikes in these data? Other thing is: its possible to identify (maybe with a dot) the pulse locations in the plot? Thanks

[Image Analyst]

How about smoothing with a Savitzky-Golay filter to get rid of oscillations then filtering out spikes that are too narrow (are noise)?

s = load('signal.mat')
signal = abs(s.signal(1:4500));
smoothedSignal = sgolayfilt(signal, 2, 101);
plot(signal);
hold on
plot(smoothedSignal, 'r-', 'LineWidth', 2);
grid on;
% Get all spikes, tall and short.
allSpikes = smoothedSignal > 0.03;
% Get rid of spikes less than 50 elements wide.
allSpikes = bwareaopen(allSpikes, 50);
[labeledSpikes, numberOfSpikes] = bwlabel(allSpikes);
props = regionprops(labeledSpikes, smoothedSignal, 'Centroid', 'MaxIntensity');
allCentroids = [props.Centroid];
xCentroids = allCentroids(1:2:end)
% Get the max intensity in each spike.
maxIntensities = [props.MaxIntensity];

% Find tall spikes
tallSpikes = maxIntensities > 0.1;
% Find short spikes
shortSpikes = ~tallSpikes;
% Plot a B at the top of every big spike
for k = 1 : length(xCentroids)
  thisX = xCentroids(k);
  thisY = smoothedSignal(round(thisX))+0.02;
  if tallSpikes(k)
    % It'a a tall spike.
    text(thisX, thisY, 'Tall', 'Color', 'r', 'FontSize', fontSize);
  else
    % It'a a short spike.
    text(thisX, thisY, 'Short', 'Color', 'r', 'FontSize', fontSize);
  end
end

meanSpacingTall = mean(diff(xCentroids(tallSpikes)))
meanSpacingShort = mean(diff(xCentroids(shortSpikes)))

[Thiago Petersen]

Hi Image Analyst, You're awesome. Your code works great for me and now i think i can continue the analysis. I'm really grateful! Cheers!