matlab实现正弦内插算法（低通滤波）

希言自然也

已于 2022-12-14 20:46:40 修改

阅读量6.4k

点赞数 5

分类专栏： MATLAB 信号处理原理文章标签： matlab 算法 python

于 2022-08-12 20:24:32 首次发布

本文链接：https://blog.csdn.net/yindq1220/article/details/126293128

版权

信号处理原理同时被 2 个专栏收录

13 篇文章 6 订阅

订阅专栏

MATLAB

6 篇文章 1 订阅

订阅专栏

一、正弦内插概述

正弦内插算法，实质是对抽样信号的还原，利用低通滤波器的时域曲线对离散信号进行插值。

根据奈奎斯特定律，抽样频率大于原始信号的最大频率的两倍以上，则可以依据抽样后信号成功还原原始信号。而数字信号处理课程中学到，时域的离散化对应频域的周期化。通过低通滤波器，截取频域上的一个周期的函数，进行离散傅里叶反变换，即可还原原始波形。

二、正弦内插公式的推导

低通滤波器在频域上的图形表示如下，横坐标Ω为其数字角频率：

时域上的函数图形表示如下：

时域上的表达公式为（傅里叶反变换）：

${\color{Red} h[t]=\frac{1}{2\pi}\int_{-\pi}^{\pi}H_{LP}(e^{jw t})dw=\frac{1}{2\pi}\int_{-w_{c}}^{w_{c}}(e^{jw t})dw =\frac{1}{2\pi}\int_{-w_{c}}^{w_{c}} (cosw t-jsinw t)dw$

对上式求积分可得：

$h[t] = \frac{1}{2\pi}\int_{-w_{c}}^{w_{c}} (cosw t-jsinwt)dw=\frac{1}{2\pi}(\frac{1}{t}sinwt+j\frac{1}{t}cosw t)\left.\begin{matrix} \\ \end{matrix}\right|\begin{matrix} \\w_{c} \\ \\ -w_{c} \\ \\ \end{matrix} =\frac{1}{2\pi t}( sinw_{c} t +jcosw_{c} t + sinw_{c} t -jcosw_{c} t) ) = \frac{sinw_{c} t}{\pi t}$

写成sinc函数形式，即为： ${\color{Red} h[t]= \frac{sinw_{c} t}{\pi t} = \frac{w_{c}}{\pi}Sa(w_{c}t)}$ 备注： $Sa(t)=\frac{sin(t)}{t}$

低通滤波器，在频域上，为原始信号的频域响应与低通滤波频域响应函数的乘积，在时域上即为原始信号与低通滤波函数h[k]的卷积。

即y[t]=x[t]*h[t]，其中x[t]为需要插值的函数，h[t]为低通滤波器的时域响应函数，y[k]为滤波后的函数，即插值后的函数。卷积公式: $y[t]=x[t]*h[t]=\int_{-\infty }^{\infty }x[\imath ]\cdot h[t-\imath]d\imath$

由于从正负无穷大内的卷积在物理上不可实现，假设插值后的离散信号的采样点的间隔周期为T，取前后合计N个离散采样点进行滤波，完成插值；假设n为离散采样点的序号，则t=nT（需要注意，x[t]的原始点保持不变，其余点置为0），因此将上述公式进行转化：

$y[t]=x[t]*h[t]=\int_{-\infty }^{\infty }x[\imath]\cdot h[t-\imath ]d\imath =\sum_{n=0}^{N-1}x[nT]\cdot h[t-nT]$

代入 $h[t]= \frac{sinw_{c} t}{\pi t} = \frac{w_{c}}{\pi}Sa(w_{c}t)$ 可得：

${\color{Red} y[t]=\sum_{n=0 }^{N-1}x[nT]\frac{w_{c}Tsin(w_{c}(t-nT))}{\pi w_{c}(t-nT)}=\sum_{n=0 }^{N-1}x[nT]\frac{sin(w_{c}(t-nT))}{\frac{\pi}{T} (t-nT)}}$

我们知道采样频率 $fs=1/T$ ，采样角频率 $w_{s}=2\pi f_{s}=2\pi/T$ ，而根据采样定理， $w_{s}>=2w_{c}$ ，因此 $wc<=w_{s}/2 = pi/T$ 。

因此有

当 $w_{c} = \pi/T$ 时，公式写成下式：

上述公式中，仅有2个参数：采样间隔T和插值点数n。（个人认为 $w_{c} < \pi/T$ 的其余值也可，此处按照 $w_{c} = \pi/T$ 简单，不同的wc即选取了不同的低通滤波器而已）

三、MATLAB正弦内插算法实现

依据此公式，使用matlab实现正弦内插的代码如下：

%sinx/x正弦内插算法实现
a=1;
%首先模拟一个离散的正弦函数
discrete_NUM = 20;%离散函数的点数共计20个
fm=zeros(discrete_NUM,1);
for m=1:1:discrete_NUM
    fm(m)=sin(0.2*pi*m);%生成正弦离散函数
end
subplot(2,1,1);
stem(fm);%正弦离散函数画图展示
hold on;

%正弦内插实现
interpolation_ratio = 10;%每个离散点插入10个点
after_NUM = discrete_NUM * interpolation_ratio;%插值后的波形的点数,一共200个点
y_o=zeros(after_NUM,1);%输出函数的幅度值矩阵
t_o=zeros(after_NUM,1);%输出函数的时间值矩阵
for t=0:1:after_NUM-1 %t取值范围为1~200
    if(mod(t,interpolation_ratio)==0) %原始离散数据位置，不计算,直接赋值
        t_o(t+1)=1+t/interpolation_ratio;%时间值，横坐标，t取值范围是1~200
        y_o(t+1)=fm(1+t/interpolation_ratio);%最终输出的幅度值 t取值范围是1~200
    else
        y=0;
        for m=1:1:discrete_NUM
            x= pi*(t-m*interpolation_ratio)/interpolation_ratio; % ，映射离散函数的采样周期为interpolation_ratio，t为插值对应的时间，每个单位映射为1
            if(x==0) %x=0时，sinx/x无意义，赋值1
                y=y+fm(m)*1;
            else
                y=y+fm(m)*sin(x)/x;%计算出每一个点内插数值
            end
%或者表达为sinc(x)的形式
%             x= (t-m*interpolation_ratio)/interpolation_ratio; % ，映射离散函数的采样周期为interpolation_num，t为插值对应的时间，每个单位映射为采样周期的1/interpolation_num
%             y=y+fm(m)*sinc(x);%计算出每一个点内插数值

        end
        t_o(t+1)=t/10;%时间值，横坐标,t取值范围是1~200
        y_o(t+1)=y;%最终输出的幅度值，纵坐标 t取值范围是1~200
    end
end
    
subplot(2,1,2);
stem(t_o,y_o);%显示插值之后的波形
hold on;

运行结果如下图。上侧为原始的离散函数，下侧为正弦内插插值后的波形。