Received signal not on sampling time

Hello all,

I have a sequence d of length N samples each of duration Ts. The channel I have will scale the received samples by a factor (1+a), which means that the received samples will not be on multiples of Ts but on multiples of Ts*(1+a). How can I capture this effect in MATLAB?

Thanks in advance

Answer 1

[Star Strider]

I’m not certain I completely understand what you want to do, but the Signal Processing Toolbox function resample may work.

[S. David]

Thanks for replying.

I don't want to resample a signal, but instead to transmit a vector of samples over a channel that scale the sample duration only. In other words, given a sequence d of N samples each of duration Ts, how to generate a received sequence of N samples each of duration Ts*(1+a)?

Thanks

[Star Strider]

I’m still not certain I understand, so see if this does what you want:

s = repmat([zeros(1,10)  ones(1,10)  zeros(1,10)],1,5);
t = linspace(1,100,size(s,2));
a = 0.25;
t2 = t*(1+a);

figure(1)
subplot(2,1,1)
plot(t,s)
axis([min(t2)  max(t2)    ylim*1.1])
title('Sent')
grid
subplot(2,1,2)
plot(t2,s)
axis([min(t2)  max(t2)    ylim*1.1])
title('Received')
grid

If it doesn’t, please explain how it doesn’t.

[S. David]

OK, let us begin again. Suppose I have the following code segment

    Ts=10^-4;
    N=512;
    d=rand(N,1);
    t=0:Ts:N*Ts;

and I want to make the sequence d aligned with sequence t, how can I do that? Can I use something like:

dt=interp1(t,d,t,'linear');

? If so, I can say with scaling it will be then:

dt=interp1(t,d,t*(1+a),'linear');

But is this true?

[Star Strider]

I don’t understand t*(1+a). To me, it looks like you simply want to expand the signal, essentially what I did.

Do you want to sample the signal (1+a) more frequently, less frequently, (both look like simple resampling to me), or shift the signal by (1+a)?

[S. David]

No, in wideband channels like underwater acoustic channels, the channel compresses/dilates the signal by the scaling factor 1+a.

[Star Strider]

I understand about frequency = propagation velocity/ wavelength, and that different water densities would affect both, but I have no idea which of these (frequency or wavelength or both) dilates or compresses.

What about this is not just resampling on a different time base to account for changing frequency or wavelength? It seems to me that the frequency would increase (in denser water) by (1+a) and the wavelength would decrease by 1/(1+a). Upsample with interpolation for lower frequencies, or downsample with decimation for higher frequencies, and you’re there. The resample function will do either of these.

[S. David]

Actually, a=v/c where v is the relative velocity between transceivers and c is the sound speed underwater. So, if the transmitted signal is s(t), the received signal then will be proportional to s(t[1+a]). So, it is more time compression/dilation rather than simple frequency shift.

Hence, I need to generate this received signal s(t[1+a]) and then do resampling at t/[1+a] to counteract the compression/dilation before detection.

[S. David]

Any help?

[Star Strider]

I’m having a very difficult time understanding what you want to do. Posting a plot of the transmitted an received signals would help a lot. I have some expertise in signal processing, but none in underwater acoustics.

[S. David]

Unfortunately, I don't have any plot of what in mind, simply because I just know the math but not the MATLAB implementation.

Are you familiar with digital communications, how we generate digital data at specific rate fs=1/Ts, convert them to analog, transmit them over a communication channel, and then sampling it at the same rate fs to get the transmitted digital data?

[Star Strider]

I’m familiar with digital communications, but haven’t done that sort of simulation in nearly 20 years. It still sounds like a resampling problem to me, especially if you know what a is.

[S. David]

There will be resampling at the receiver on the received signal. Do you mean we can generate the received signal itself using resampling?

[Star Strider]

If I understand what you’re doing (no guarantees), it seems to me that resampling would regenerate your transmitted signal, or with the appropriate dilation or compression, simulate your received signal. That was my first understanding of your description of the process you’re analysing.

[S. David]

How to do this in MATLAB given a transmitted vector d and a time vector t?

[Star Strider]

If you have the Signal Processing Toolbox, use the resample function. There are also links to the documentation of a number of related functions at the end of that page. (I hyperlinked to the online documentation for it here.)

[S. David]

Thanks

[Star Strider]

Did it solve your problem?

[S. David]

I saw the resample function in MATLAB, and it says that it " resamples the sequence in vector X at P/Q times the original sample rate using a polyphase implementation ". But the question is how set P and Q to resample my signal at t*(1+a)? I think P/Q must equal 1+a, but a is not an integer!!

[Star Strider]

Here’s an example using interp1. I’m still not certain I understand what you’re doing, but this may help. Please experiment with it. Does it do what you want?

t = linspace(0, 4*pi, 35);
s = sin(t)+cos(3*t);

figure(1)
plot(t, s)
grid

t2 = linspace(0,5*pi,50);
s2 = interp1(t, s, t2, 'pchip', 'extrap');

figure(2)
plot(t2, s2)
axis([xlim    -2  2])
grid

[S. David]

I think it will resemble my case more if the time range changes while the number of points be fixed. Like this:

t = linspace(0, 4*pi, 35);
s = sin(t)+cos(3*t);

figure(1)
plot(t, s)
grid

t2 = linspace(0,4*pi*(1+0.024),35);
s2 = interp1(t, s, t2, 'pchip', 'extrap');

figure(2)
plot(t2, s2)
axis([xlim    -2  2])
grid

I think your example in your second reply is clearer since it plots rectangular pulses. But there, you didn't align s to t2 as here.

[Star Strider]

OK, change the first linspace length to 36 and :

s = [0 0 repmat([zeros(1,4) ones(1,4)], 1, 4) 0 0];

This is a square-wave sequence for you to experiment with. Change the interp1 method to 'linear' if you want. The method is probably not critical with this signal.

I’m simply trying things randomly since I still don’t understand what you’re doing. Everything I’ve suggested so far has gone for naught.

[S. David]

Basically I want to transmit an K-by-1 vector d over a channel with linear time varying impulse response:

h(tau;t)=sum_p h_p delta(tau-[tau_p-a*t])

The passband transmit signal is given by:

s(t)=Re{exp(j2*pi*fc*t) sum_k d_k g(t-kTs)}

where fc is the carrier frequency and g(t) is a rectangular pulse of duration Ts (the symbol time).

The baseband received signal in this case is given by:

v(t)=sum_k d_k sum_p h_p exp(-j2*pi*fc*tau_p) exp(j2*pi*fc*t*a) g(t[1+a]-kTs-tau_p)+w(t)

where w(t) is additive white Gaussian noise. So, you can see that the symbol duration has changed from Ts to Ts/(1+a) from the argument of g(t) in the transmit and receive signals. Eventually, I want to generate the signal v(t) from the vector d.

That is what I am trying to do. I do it in MATLAB like this:

for pp=1:P
    yT=interp1(t,d,t*(1+a)-tau(pp),'cubic');
    y=y+yT*h(pp)*exp(-1i*2*pi*fc*tau(pp)).*exp(1i*2*pi*fc.*t*a);
end

Is this correct?

[Star Strider]

It is if it runs and produces the result you want.

[S. David]

Actually, what I do is that I compare y above with its numerical equivalent, and they don't match, where mathematically y=Hx, and I know H and x!! But there is a chance that forming H has a problem!!

[Star Strider]

Your h looks like a convolution, and s(t) like a cosine transform of some sort?

By the way, resample simply states that p and q must both be positive integers. You may be able to use the core MATLAB rats function to create a ratio of two integers for 1+a or 1/(1+a) that resample will accept.

To get the numerator and denominator from the rats function to use in resample, this will work:

x = 1.15077
y = rats(x);
xnumden = [y(1:findstr(y,'/')-1) ' '  y(findstr(y,'/')+1:end)];
xnd = str2num(xnumden)

[S. David]

Interesting!! So, I can write

x=(1/(1+a));
.
.
.  
d_res=resample(d,xnd(1),xnd(2))

to get the time-scaled version of d, right? But now what about the delay shift effect tau? How to capture it in this way?

[Star Strider]

That would do for the time-scaling. I got the impression that tau was involved in the convolution. Not sure how you want to deal with that.

[S. David]

OK, I will figure it out. Thanks for your patience and help.

[Star Strider]

My impression is that tau was a variable of integration in a convolution, so it never appeared in the integrated (convolved) result, only time t. If you’re simply shifting it in time by tau, first resample, then shift.

[S. David]

If we ignore the time scaling, the noise-free baseband received signal would be:

v(t)=sum_k d_k sum_p h_p exp(-j*2*pi*fc*tau_p) g(t-kTs-tau_p)

sampling this at t=nTs results in:

v(nTs)=sum_k d_k f_{n-k}

where

f_m=sum_p h_p exp(-j*2*pi*fc*tau_p) g(mTs-tau_p)

for m=0, ..., L, and L is the maximum delay in symbol time.

So, the discrete time equivalent signal can be written as:

v_n=sum_{m=0}^L f_m d_{n-m}

which is the convolution.

However, if I incorporate the time scaling, then

v_n=v(nTs)=sum_k d_k f_{n,n-k}

where

f_m=sum_p h_p exp(-j*2*pi*fc*tau_p) exp(j*2*pi*fc*n*Ts*a) g(mTs-[tau_p-a*n*Ts])

So, it appears to be convolution, however, each n has its own channel length. In this case, how can I make convolution for each received symbol separately? Is there any built-in function to do this?

[Star Strider]

There are convolution functions in core MATLAB and the Signal Processing Toolbox, that also has modulation and demodulation functions. I don’t know if they’re relevant or not.

There might also be in the Communications System Toolbox. I don’t have it, so I haven’t explored it. (I do mostly simulation, modeling, parameter estimation, control, and signal processing.)

[S. David]

OK, I will check.

You mentioned that I can shift the resampled signal in time. How?

[Star Strider]

To shift it in time as I was thinking of it means generating two zeros vectors both of length t+dt, where t is the signal length and dt the length of the time shift. Copy the same signal to one vector starting at t=0 and to the second vector starting at t=t+dt, so the first one starts at t=0 with a gap of dt at the end, and the second has a gap of dt at the beginning and goes to the end. (It creates an echo effect with audio files.)

[S. David]

I see. But in this case I need to convert the tau (ms) to symbol time, right?

[Star Strider]

I assume so. I was describing it in terms of vector indices because that’s how I would code it.

[S. David]

Actually, I see what you were saying now. So, I don't need the time vector in this case. I define the transmit vector as:

dCP=[d(N-L+1:end) d];%Cyclic prefixed signal

Then I find the resampled signal as following:

a=0.004;

x = (1+a);
s = rats(x);
xnumden = [s(1:findstr(s,'/')-1) ' '  s(findstr(s,'/')+1:end)];
xnd = str2num(xnumden);
p=xnd(1);
q=xnd(2);

v=resample(dCP,p,q)

Now for the time shift, I do the following:

y1=[v zeros(1,L)];
y2=[zeros(1,L) v];
y=y1+y2;

where I assumed I have two time shift one at tau=0 and one at tau=L (maximum delay in symbol time). After finding the time shift, I do resampling again to recover the original data as:

x = 1/(1+a);
s = rats(x);
xnumden = [s(1:findstr(s,'/')-1) ' '  s(findstr(s,'/')+1:end)];
xnd = str2num(xnumden);
p=xnd(1);
q=xnd(2);

yRes=resample(y,p,q)

Do you think I am interpreting it right?

[Star Strider]

It looks correct to me. If it produces results that make sense in the context you’re applying it, then we have likely arrived at a solution!

[S. David]

I need to check it. But in this way, where is the effect of Ts (symbol time) and Ta (sample time)? I mean I have a bandwidth B, Ts=1/(2*B), and Ta=Ts/1000.

[Star Strider]

We’re beyond my expertise in digital communications. I very much hope it works!

[S. David]

You are awesome!! Using this approach my code works very well. I tried two cases: 1-) when a=0 (there is no time scaling), and 2-) when a=0.004. In both cases the noise-free received signal exactly matches its numerical equivalent.

However, when a~=0, the received signal contains the factor exp(1i*2*pi*fc*a*t). Now since I am not using t, how can I incorporate it in the received signal term?

Note: In the receiver, the received signal is multiplied by exp(-1i*2*pi*fc*a*t), to counteract the effect of Doppler shift, so my code works fine without incorporating exp(1i*2*pi*fc*a*t) in the received signal. However, I want to know how, because I will need it for more complicated channel models!!

[Star Strider]

Thank you!

The ‘extra’ factor seems to me a phase term of the sort I would expect as phase delay in analog filter output (other than a Gauss filter). I’m casting it in terms of my experience, so filters come quickly to mind. Does the received signal have characteristics of a filtered signal?

(Which brings to mind the possibility that you could consider your channel as a filter?)

[S. David]

The channel is a Linear Time Varying filter with an impulse response

h(tau;t)=sum_p h_p delta(tau-[tau_p-a*t])

The passband transmitted signal is given by:

s(t)=Re{exp(j*2*pi*fc*t) sum_k d_k g(t-kTs)}

Passing s(t) through h(tau;t) (the convolution between both of them) and downconvert the signal to baseband with respect to the carrier frequency fc, the noise-free baseband received signal is given by:

v(t)=sum_k d_k sum_p h_p exp(-j*2*pi*fc*tau_p) * exp(j*2*pi*fc*t*a) g(t(1+a)-kTs-tau_p)

That is how the factor exp(j*2*pi*fc*t*a) is appeared!!

[Star Strider]

I’m doing my best to follow this. So demodulating it introduced the phase shift? That would seem reasonable.

[S. David]

Passing the signal through the channel causes the phase shift. I've done it this way: First I defined the time vector t as

t=CP*Ts:Ts:(N+CP-1)*Ts

and then after time scaling and shifting the received signal y is multiplied by exp(1i*2*pi*fc*a.*t).

Thank you very much. This discussion was of a great help, I really appreciate it, and I hope we can continue discussion as I go, since although I am good at math, I am not that good in interpreting math to a MATLAB code.

[Star Strider]

My pleasure!

It was an education for me, too. My last communications course was a bit over 15 years ago (my interest being biomedical telemetry), and I haven’t done much with it since. I had to search through my library to find my textbooks and review them. I vaguely remember homework problems on the situation you’re dealing with, but obviously not in any great detail. I did discover my ancient MATLAB code for those assignments.

I’ll always be happy to provide whatever expertise I may have.

[S. David]

I think it is not necessary to be familiar exactly with the material. I am sure an accurate description of what I want to do will suffice for you to suggest how to do it.

Thanks again!!

[Star Strider]

My pleasure!

I’m happy to help.

[S. David]

Hi again,

When I have one factor a, the channel multiplies the signal by exp(j*2*pi*fc*t*a) and the receiver counteracts this by multiplying by exp(-j*2*pi*fc*t*a). So, I didn't have to include this in my code, since mathematically the effect of the channel and receiver processing cancels each others.

Now I am trying to resample the signal several times, where the received signal now is written as:

v(t)=sum_k d_k sum_p h_p exp(-j*2*pi*fc*tau_p) * exp(j*2*pi*fc*t*a_p) g(t*(1+a_p)-kTs-tau_p)

How can I deal with this in MATLAB?

t is defined as

t=0:Ts:(N-1)*Ts;

which has the same length as d. However, after resampling d will have p/q times the length of d!! I tried to do the following:

for cc=1:Np
    x1 =(1+a(cc));
    s1 = rats(x1);
    xnumden1 = [s1(1:findstr(s1,'/')-1) ' '  s1(findstr(s1,'/')+1:end)];
    xnd1 = str2num(xnumden1);
    p1(cc)=xnd1(1);
    q1(cc)=xnd1(2);
end

Ld=max(p1./q1)*N;

y=zeros(1,Ld);

for cc=1:Np
z=resample(d,p1(cc),q1(cc));
r=p1(cc)/q1(cc);
TsNew=Ts/r;
tNew=0:TsNew:(N-1)*Ts
y(1:length(z))=y(1:length(z))+z.*exp(1i*2*pi*fc*a(cc).*tNew);
end

The problem is tNew and z's sizes always differ by one!!

1- Is this approach correct in principle?

2- Why tNew and z are not of the same size as they should theoretically?

Thanks

[Star Strider]

This is not my area of expertise, so I can only look at your code to see if I can find anything that might be a problem.

It will probably be best for you to define tNew by its beginning and end values (that you calculate) and then use the linspace function to define it to have the same length as z. That may be the most efficient (and least agonizing) way to equalise their lengths. I doubt if it will lead to significant loss of accuracy or precision in your subsequent signal processing.

[S. David]

You mean like this:

tNew=linspace(0,(N-1)*Ts,length(z));

right? I guess it makes sense. I will continue using this and see what happens.

Thanks

[Star Strider]

Yes. It may be ‘cheating’ a bit, but sometimes that’s necessary.

[S. David]

Actually, the step size in tNew (tNew(2)-tNew(1)) is the same as Ts/r where r=length(z)/N. I think it is accurate enough.

[Star Strider]

If z is relatively long, there may not be much difference from lengthening or shortening tNew by one element. The key is the result: the lengths are equal.

[S. David]

I agree.

Answer 2

[rifat]

I think the following script should do the work:

signal_vector= [1 1 1 1 1]; % Define your signal here
Ts=1/50; % Define sampling period here
a=1; % define 'a' here
N=length(signal_vector);
x=0:Ts:(N-1)*Ts;
new_x=0:(1+a)*Ts:(N-1)*(1+a)*Ts;
new_signal=interp1(x,signal_vector,new_x,'cubic'); % change the algorithm flag to 
                                                   % 'pchip' or 'spline' if
                                                   % you need

[S. David]

Thanks for replying.

I wrote something like this in the third post to the first answer, and I see it is basically the same, except the parameter 'cubic', I had it 'linear' in the interp1 function. Does it matter?

[rifat]

Sorry.. I didnt read other answers and comments...

[S. David]

That is OK, but is there any significant difference between the flags of the interp1 function, or all serve the same end in different ways?

[Star Strider]

There are. The method you use depends on what you’re interpolating and how you want to do it.

[rifat]

That depends on how much smoothness you want in the interpolated signal. For less sophisticated applications, from my experience, the results dont vary much in different flags. But, avoid using 'linear' unless you want a very coarse interpolation.

[rifat]

Seems like you are working some kind sonar application. I worked on ultrasound imaging that has almost same philosophy and also needed same kind of stretching/compression operation to account for tissue motion. I used the program above.

[S. David]

For example, for the problem in hand, why did you use 'cubic' and not 'linear'?

[rifat]

If you use linear, then actually the samples are connected by a straight line and intermediate values are obtained from the fitted straight line. In other cases, some sophisticated smooth curves are fitted to obtain the intermediate values. So, 'linear' will, in general, result in more noise.

[S. David]

Thanks

[John D'Errico]

Note: I once did an analysis where I showed that IF you have noise in the process, so you are interpolating a combination of noise PLUS signal, then linear is often a better answer, depending on the signal to noise ratio. Cubic interpolation on the noise component has a higher variance than does linear.

[S. David]

No, I am adding the noise after interpolation. So, in this case, 'cubic' works fine I guess. But I am not sure how accurate the whole interpolation process.