how do I re-sample a continuous function; but maintain shape?

10 次查看（过去 30 天）
Eric Kuebler
Eric Kuebler 2019-11-5
评论： Eric Kuebler 2019-11-6
I am trying to resample a vector that represents a voltage time-series from 200kHz to 50 kHz using the resample function (the two vectors are shown below).
untitled.jpg
The code I've used to resample is as follows:
x = resample(x,5e4,20e4);
When I measure what we call the 'half width' (the duration of time at half the height of the voltage peak) I notice that there is an increase, essentially the half width becomes longer. The increase also varies between waveforms.
I am curious if there is anyway to reduce this discrepancy, or at least perhaps have it consistent and based on the slower sampling rate? I've tried using filter coefficients that mimick the shape of these waveforms; but that seemed to make things worse.
Thank you in advance!
Eric
  5 个评论
显示 3更早的评论
Eric Kuebler
Eric Kuebler 2019-11-5
编辑：Eric Kuebler 2019-11-5
My apologies, I tried downsampling to 10kHz and this is where the error is much much larger.
Comparing the waveforms above (so 200kHz to 50kHz) there is a 7% difference.
Eric Kuebler
Eric Kuebler 2019-11-5
I am computing the half width by the following
dt = 1000/5e4;
[peakV,peakTime] = max(x); % peak
[troughV,troughTime] = min(x(peakTime:end)); % trough
troughTime = peakTime+troughTime-1; % adjusted trough time
halfHeight = (peakV-troughV)/2; % half height
halfHeight = peakV-halfHeight; % adjusted half height
halfHeightTpre = find(x(1:peakTime)<halfHeight,1,'last'); % half height pre-peak
halfHeightTpost = peakTime+find(x(peakTime:troughTime)<halfHeight,1,'first'); % half height post-peak
hw = (halfHeightTpost-halfHeightTpre)*dt; % half-width

请先登录，再进行评论。

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

采纳的回答

Daniel M
Daniel M 2019-11-5
编辑：Daniel M 2019-11-5
I see several potential reasons why this is occuring.
1) There will be general broadening of the signal due to downsampling.
2) Since the magnitude of the downsampled signal is potentially lower, then half of the max could also be lower. So you could be measuring a wider part of the signal.
3) You're measuring half-width from the left side. This combined with issue 2 could lead to measuring something like the blue line in "Capture" (attachment). I would have defined the half-width as something more like the red line, which is the HWHM from the left side.
4) You're using discrete values to calculate width. Probably not an issue because the sampling rates are still quite high, but could lead to something occuring like in "Capture2" attachment. If you find this is a big issue, then you could interpolate between the first point <halfHeight and first point >halfHeight.
I'm not really sure what to suggest other than to try using the functions pulsewidth and/or statelevels in the signal processing toolbox.
  1 个评论
Eric Kuebler
Eric Kuebler 2019-11-6
  1. I agree, this is why I mentioned that perhaps there is a way to do this that the difference would be based on the slower sampling rate.
  2. I agree, if this is the case I would expect the same amount of error in a trace that was originally sampled at 50kHz.
  3. It may be useful for us to try taking this measure from the left side (i.e., threshold).
  4. This could also be the case; but agree sampling rates are still quite high.I will assess whether this is an issue.
Can you say a bit more about how I might use the function statelevels?
BTW, I didn't mean to accept this answer, I was hoping to leave it open; but hit that button by accident, good thing they put the 'Accept' button so close to the attachments!!

请先登录，再进行评论。

更多回答（1 个）

Daniel M
Daniel M 2019-11-6
You can unaccept the answer if you want to keep it open, but also consider this answer.
(For 3, I meant you're measuring from right side, and could try from left. I think you understood what I meant).
Here are some ideas. I'm not saying they are robust, but this is certainly a potential direction to go in.
% create sample data
fs = 4e6;
t = (0:99).*1/fs;
y = [ones(1,45) 1+sin(2*pi*320000*t(1:10)) zeros(1,45)+2e4.*t(1:45) ] + 0.1.*rand(1,length(t));
% that looks ridiculous I know..., but it isn't a bad approximation to your signal
% get the statelevels
levs = statelevels(y);
% if you call it without the output, it will display the attached figure "levels"
% get the overshoot value, need to flip y
% I'm sure using max or findpeaks is sufficient here too.
% I'm just showing you the utility of these functions
[os,oslev,osinst] = overshoot(fliplr(y),t,'StateLevels',levs);
% if you plot without outputs, get attached figure "overshoot"
% we care about 'upper state' and 'post-overshoot'
% get the pulse width, the crossings, and the midlevel
[width,init,fin,midl] = pulsewidth(y,t,'StateLevels',[levs(2) oslev]);
% width = fin-init
% see attached figure "pulsewidth"
As you can see, these functions can be pretty useful! I believe these functions have the ability to work on trains of pulses too, so you don't even need to epoch your data.
If nothing in here quite works for you, then I am curious what the application of your problem is for. Maybe there is another way to circumvent the issue.
  1 个评论
Eric Kuebler
Eric Kuebler 2019-11-6
I am going to try this in our analysis, and I appreciate the time you 've taken to make the advantages of these functions clear to me.
It is important for our application that the function choses either pre- or post-peak as the lower state consistently. I think that is the only advantage that the code I loaded above may have, I tell it to look in the post-peak portion for that minimum.
The application we are using this for is classification of cell types in single unit recordings of neurons in the visual cortex. The half width duration is one of the parameters that gives good classification.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Multirate Signal Processing 的更多信息

标签

产品


版本

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by