Code for signal analysis

15 次查看（过去 30 天）
Thiago Petersen
Thiago Petersen 2016-8-13
评论： Thiago Petersen 2016-8-14
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);
  2 个评论
John BG
John BG 2016-8-14
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:jgb2012@sky.com jgb2012@sky.com>
Thiago Petersen
Thiago Petersen 2016-8-14
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!

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Image Analyst
Image Analyst 2016-8-13
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)
  3 个评论
显示 1更早的评论
Image Analyst
Image Analyst 2016-8-14
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
Thiago Petersen 2016-8-14
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!

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Spectral Measurements 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by