采样率转换理论基础

最新推荐文章于 2023-07-19 22:43:45 发布

半截木头渡海洋

最新推荐文章于 2023-07-19 22:43:45 发布

阅读量1w

点赞数 11

分类专栏： DSP 好玩的东西数学基础文章标签：插值重采样

本文链接：https://blog.csdn.net/ljl86400/article/details/81566879

版权

DSP 同时被 3 个专栏收录

19 篇文章 2 订阅

订阅专栏

数学基础

16 篇文章 6 订阅

订阅专栏

好玩的东西

3 篇文章 0 订阅

订阅专栏

离散信号重采样（采样率转换）通常有两种实现办法

1】一种是通过DA转换成模拟信号然后再通过AD重新对模拟信号进行采样

2】一种是通过插值的办法在对输入信号进行插值或者抽取的方式，完成纯数字域的采样率转换

本文主要记录第二种方法

1.通过插值公式实现对原始信号的重构

$g(t)=\frac{\sin (\pi t/ T)}{\pi t/ T}$

$y(t)=\sum _{n=-\infty}^{\infty}x(nT)g(t-nT)$

观察以上公式，gt函数中sin(pi*t/T)的周期为1/2*Fx=1/(2T)，即gt频率是原信号采样频率的1/2。

< 例一 >

先形象的了解一下插值函数g(t)

假设我们有信号连续时间信号x（t）=cos(2*pi*f*t+phi)，并以Fs=20的速度对x(t)进行采样，得到x(n)，并有如下图形：

现在我们要根据采样获得的这些离散点，结合插值的方法得到原来的xt连续信号：

在上图中，经过n=0:20的累加过程，已经有一部分信号十分贴近于连续原始信号xt，只是在两端点位置因为截断问题，有一些频谱泄露，可以证明如果n的取值范围为(-∞，∞)时，可以完整重构原始连续信号xt。

2.离散时间信号的插值重构

现在我们对插值公式进行一些处理，对插值信号gt按照1/Ty的频率进行采样，就可以得到一个重构函数序列的计算公式

$g(mT_{y})=\frac{\sin (m\pi T_{y}/T)}{m\pi T_{y}/T} \qquad m=0,1,2,3...$

序列结构如下图

将采样周期为Tx的输入x（nTx）序列与这个插值序列卷积，就得到了一个经过插值的序列y（mTy），插值后采样周期变为Ty，插值后输出序列的第m个元素记为以下形式：

$y(mT_{y})=\sum_{n=-\infty }^{\infty}x(nT_{x})g(mT_{y}-nT_{x})\qquad m=0,1,2,3,...$ 公式1

如果输入序列x长度n=10,那么每计算一个y(mTy)值，需要相应的将g(mTy)计算出来后平移n-1次，然后乘以系数x(nTx)，再求加和，也就是说x序列越长，g就需要求取相应更多次的平移序列，这需要很大的计算量。

通过对公式1的观察，我们可以看出，插值后的结果y(mTy)是x(nTx)与g(nTx)卷积的结果（也可以将y(mTy)看做是线性时不变系统g(nTx)对信号x(nTx)的响应），即Y(z)=X(z)*G(z),卷积后结果y(mTy)是对原是信号xt的周期采样,也即y(mTy)=x(nTy)，采样周期为Ty，采样频率为1/Ty，由此完成了x(nTx)到x(mTy)的采样频率转换。同时我们观察会发现，插值信号的采样频率1/Ty最终决定了y(mTy)的采样周期，因此要想从1/Ty中恢复xt，1/Ty须大于等于xt中最高频率成分的2倍频。

针对公式1，进行代码验证，代码如下

clc
clear 
close all

%% 参数配置
t = -1:0.001:1;
tx1 = 0:0.001:1;

% 原始信号频率以及采样参数，1/Tx=Fx≥2fsignal,Fx≥2*fignal,Fx=2*fg
fsignal = 2;        % 原连续信号的频率,通常插值函数的频率≥2*原始信号的频率
Fx = 20;            % 对原信号进行采样的采样频率，应该≥2fsignal,对应gt的频率1/2
phi = 0.2;
Tx = 1/Fx;
tx = 0:Tx:1;


% 插值函数采样相关参数，插值函数的原始频率fg=1/2*Fx≥fh
Fsy = 5;           % 对插值函数进行采样的频率，应该大于等于2*fsignal
Ty = 1/Fsy;
txy = 0:Ty:1;
ty = (-1:Ty:1)+10^(-10);

%% 插值函数原始图形
gt = (sin(pi*t/Tx)./(pi*t/Tx));
figure(1)
subplot(311)
plot(t,sin(pi*t/Tx))
title('sin(pi*t/Tx)')
xlabel('t=-1:0.001:1')
subplot(312)
plot(t,(pi*t/Tx))
title('pi*t/Tx')
xlabel('t=-1:0.001:1')
subplot(313)
plot(t,gt)
title('gt')
xlabel('t=-1:0.001:1')
ylabel('gt')

%% 信号函数与采样
xt = cos(2*pi*fsignal*tx + phi);
xt1 = cos(2*pi*fsignal*tx1 + phi);
figure(2)
hold on
stem(tx,xt,'LineStyle','none')
plot(tx1,xt1)
title('以Fs为采样周期对xt进行采样')
xlabel('tx=0:1/Fx:1,tx1=0:0.001:1')

%% 信号重构
ytn = zeros(length(tx),length(t));
xtn = zeros(length(tx),length(t));
for iloop1 = 1:Fx+1
    ytn(iloop1,:) = xt(iloop1)*(sin(pi*(t-(iloop1-1)*Tx)/Tx)./...
        (pi*(t-(iloop1-1)*Tx)/Tx));
    xtn(iloop1,:) = t+(iloop1-1)*Tx;
end

figure(3)
subplot(411)
hold on
plot(t,ytn(1,:))
plot(t,ytn(2,:))
plot(t,ytn(3,:))
hold off
legend('n=0','n=1','n=3')
title('x(nTx)*g(t-nTx)')
xlabel('t=-1:0.001:1')
subplot(412)
hold on
plot(t,sum(ytn(1:3,:)))
plot(t,ytn(4,:))
hold off
legend('sum(1:3)','n=4')
title('x(nTx)*g(t-nTx)')
xlabel('t=-1:0.001:1')
subplot(413)
hold on
plot(t,sum(ytn(1:4,:)))
plot(t,ytn(5,:))
hold off
legend('sum(1:4)','n=5')
title('x(nTx)*g(t-nTx)')
xlabel('t=-1:0.001:1')
subplot(414)
hold on
plot(t,sum(ytn))
plot(t,ytn(Fx+1,:))
plot(tx1,xt1)
hold off
legend('sum(1:21)','n=21','xt')
title('重构信号ytvs原始连续信号xt对比')
xlabel('t=-1:0.001:1')

%% 对原连续信号xt按照周期Ty进行采样，得到xty
xty = cos(2*pi*fsignal*txy + phi);
xt1 = cos(2*pi*fsignal*tx1 + phi);
figure(4)
hold on
stem(txy,xty,'LineStyle','none')
plot(tx1,xt1)
title('以Fs为采样周期对xt进行采样')
xlabel('txy=0:Ty:1,tx1=0:0.001:1')


%% 对重构信号进行采样，采样周期为Ty
gty = (sin(pi*ty/Tx)./(pi*ty/Tx));
figure(5)
hold on 
stem(ty,gty,'LineStyle','none')
plot(t,gt)
hold off
title('对gt进行采样')
xlabel('t=-1:0.001:1,ty=-1:1/Fsy:1')


%% 利用两次采样的结果xty，gty得到一个新的信号，验证这个信号是否是原连续信
% 号按照Ty采样的结果
ytnd = zeros(length(tx),length(ty));
xtnd = zeros(length(tx),length(ty));
ytndSum = zeros(length(tx),length(ty)*(length(tx)-1));
xtndSum = zeros(length(tx),length(ty)*(length(tx)-1));
tySum = -1:Ty*Tx:1+(Fx-1)*Ty*Tx;
for iloop1 = 1:Fx+1
    ytnd(iloop1,:) = xt(iloop1)*(sin(pi*(ty-(iloop1-1)*Tx)/Tx)./...
        (pi*(ty-(iloop1-1)*Tx)/Tx));
    xtnd(iloop1,:) = ty+(iloop1-1)*Tx;
end
figure()
subplot(311)
grid on
hold on
stem(xtnd(1,:),ytnd(1,:))
stem(xtnd(2,:),ytnd(2,:))
stem(xtnd(3,:),ytnd(3,:))
stem(xtnd(4,:),ytnd(4,:))
xlim([-1 1])
title('x(nTx)*g(mTy-nTx),n=0,1,2,3')
legend('n=0','n=1','n=2','n=3')
hold off
subplot(312)
grid on
hold on
stem(ty,sum(ytnd(1:2,:)))
stem(xtnd(3,:),ytnd(3,:))
xlim([-1 1])
legend('sum(1:2)','n=2')
title('x(nTx)*g(mTy-nTx),m=0')
hold off
subplot(313)
grid on
hold on
stem(ty,sum(ytnd(1:Fx+1,:)))
plot(tx1,xt1)
xlim([-1 1])
legend('sum')
title('sum(x(nTx)*g(mTy-nTx))')
hold off

3.频率转换公式的进一步变形与探讨

假设Tx=Ty=T，我们就可以将以上公式1简化，得到一个输入信号与线性时不变系统的卷积公式

$g(m)=\frac{\sin (m\pi )}{m\pi } \qquad m=0,1,2,3...$

$y(m)=\sum_{n=-\infty }^{\infty}x(n)g(m-n)\qquad m=0,1,2,3,...$ 公式2

如果Tx≠Ty，将公式1整理得到如下形式

$y(mT_{y})=\sum_{n=-\infty}^{\infty}x(nT_{x})g(T_{x} (\frac{mT_{y}} {T_{x}}-n))$ 公式3

我们可以将mTy/Tx的结果分成整数km与小数△m两部分，即：

$\frac{mT_{y}}{T_{x}}=k_{m}+\Delta_{m}$ 公式4

将公式4带入到公式3中，得到以下形式：

$y(mT_{y})=\sum_{n=-\infty}^{\infty}x(nT_{x})g(T_{x}(k_{m}+\Delta_{m}-n))$ 公式5

如果我们定义 $k=k_{m}-n$ ，公式5可以整理为以下形式

$y(mT_{y})=\sum_{k=-\infty}^{\infty}x((k_{m}-k)T_{x})g(kT_{x}+\Delta_{m}T_{x})$

同时，mTy=（km+△m）Tx，因此：

$y(mT_{y})=y((k_{m}+\Delta_{m})T_{x})$

也即：

$\small y(mT_{y})=y((k_{m}+\Delta_{m})T_{x})=\sum_{k=-\infty}^{\infty}g(kT_{x}+\Delta_{m}T_{x})x((k_{m}-k)T_{x})$ 公式6

综合以上讨论：

1】给定的Tx决定了输入序列的采样率，给定Tx就确定了输入序列的时间坐标；输出的采样率由Ty确定，给定Ty就确定了输出序列的时间坐标

2】讨论观察g(kTx+△mTx)序列后我们能发现，g(kTx+△mTx)是g(kTx)的尺度缩放，而△m的取值可能是有限个的，也就是说在公式6中用到的g(kTx+△mTx)序列只有有限种可能，这样相比较公式1中需要求取每一个n对应的g序列来讲，我们就可以用提前计算所有可能的g(kTx+△mTx)然后再查表的方法来代替原来公式1中每次都需要重新计算序列的方法，达到降低计算量的目的；而x((km-k)Tx)部分则是原公式1中x(nTx),可以通过km与n的对应关系查表获得，无需重复计算。

3】通过以上讨论我们也能发现，对于不同的输入样本点，系统的响应函数是在变化的，所以离散时间信号的速率转换系统一个线性时变系统。

以上就是信号采样率的转换理论基础，上采样与下采样的讨论结合子带滤波器组，后续再开一篇博客继续讨论。