手撕重采样,考虑C的实现方式

news2024/11/16 15:29:06

一、参考文章:

重采样、上采样、下采样 - 知乎 (zhihu.com)

先直接给结论,正常重采样过程如下:

1、对于原采样率fs,需要重采样到fs1,一般fs和fs1都是整数哈,则先找fs和fs1的最小公倍数,设为m,设m/fs=M,m/fs1 = L。则信号先要做M倍的插值,即上采样,再做1/L倍的抽取,即下采样;

2、因为插值和抽取,信号的频带都会变,也就是信号会引入噪声,所以需要滤波处理;

3、具体来说,插值,频谱变窄,即信号频带压缩了,如果不做处理,信号会包含带宽以外的噪声,所以需要做低通只滤出变窄的信号频带,去掉噪声。抽取,之后的频谱会变宽,即最终和原信号频谱对不上,所以又需要低通滤波,将信号滤到和原信号一样的频带。

4、具体来说,第一个低通的通带(归一化)是0~1/M,第二个低通的通带是0~1/L。找到满足性能的滤波器,还是比较有考验的。

回到程序,上文中第一段,涉及到滤波器,改用0相位滤波,即调用matlab的filtfilt。滤波时考虑使用原程序中的滤波器幅值归一化。且,加入抽取、滤波,即实现重采样完整流程,最后,还对比matlab自己的resample函数结果。如下:

clc;
close all;
clear all;
%%%%%%%%%%%%%%%%%程序说明
% 1、使用'zero stuffing'和'low pass filter'实现内插/上采样

N = 4; % 上采样率
h_t = ones(1, N);
h_t_lp = sinc(-pi:2*pi/20:pi);
% h_t_lp = h_t_lp ./ sum(h_t_lp); % 归一化LPF


A = 1.5;
B = 1;
f1 = 50;
f2 = 100;
Fs = 1000;
t = 0:1/Fs:1;
sig = A * cos(2 * pi *f1 * t) + B * sin(2 * pi *f2 * t); % 原始数据

% 插值:插入0
sig_zerostuff = zeros(1, length(sig) * N);
sig_zerostuff(1 : N : length(sig_zerostuff)) = sig;
sig_zerostuff_lp = conv(sig_zerostuff, h_t_lp);
% 采样保持:使得插入的数值与原数据的数值一致
sig_sample_hold = conv(sig_zerostuff, h_t);
% 将采样保持的信号经过LPF
h_t_lp = h_t_lp ./ sum(h_t_lp); % 归一化LPF
% sig_sample_hold_lp = conv(sig_sample_hold, h_t_lp);
sig_sample_hold_lp = filtfilt(h_t_lp,1,sig_sample_hold);

subplot(5,1,1);
stem(sig,'MarkerFaceColor',[0 0 1]);xlim([1 25]);
title('原始数据');

subplot(5,1,2);
stem(sig_zerostuff,'MarkerFaceColor',[1 0 0]);xlim([1 25*N]);
title('插值后的数据');

subplot(5,1,3);
stem(sig_zerostuff_lp,'MarkerFaceColor',[1 0 0]);xlim([1 25*N]);
title('插值后的数据通过LPF');

subplot(5,1,4);
stem(sig_sample_hold,'MarkerFaceColor',[1 0 0]);xlim([1 25*N]);
title('插值保持后的数据');

subplot(5,1,5);
stem(sig_sample_hold_lp,'MarkerFaceColor',[1 0 0]);xlim([1 25*N]);
title('采样保持的数据通过LPF');

% 抽取,1/3
sig_sample_hold_lp_downsample = sig_sample_hold_lp(1:3:end);
% sig_sample_hold_lp_downsample_lp = conv(sig_sample_hold_lp_downsample, h_t_lp);
sig_sample_hold_lp_downsample_lp = filtfilt(h_t_lp,1,sig_sample_hold_lp_downsample);
sig_resample = resample(sig, 4, 3);

figure
subplot(4,1,1);
stem(sig,'MarkerFaceColor',[0 0 1]);xlim([1 120]);
title('原始数据');
subplot(4,1,2);
stem(sig_sample_hold_lp_downsample,'MarkerFaceColor',[0 1 0]);xlim([1 120]);
title('抽取数据');
subplot(4,1,3);
stem(sig_sample_hold_lp_downsample_lp,'MarkerFaceColor',[1 0 0]);xlim([1 120]);
title('抽取数据通过LPF');
subplot(4,1,4);
stem(sig_resample,'MarkerFaceColor',[0 0 0]);xlim([1 120]);
title('resample重采样数据');

wAAACH5BAEKAAAALAAAAAABAAEAAAICRAEAOw==

结果:

还是文章一,插值可以用另一种方法,即,不是插0值或者保持,而是用其他插值法,matlab有interp函数,可以设置参数用不同插值法。典型的线性插值,应该比较好找C代码。

先看看文章的效果:

效果还是可以的。

关于线性插值C代码,时间关系暂未验证哈:

线性插值_c语言实现_线性插值c语言程序-CSDN博客

标准C语言插值函数 - 百度文库 (baidu.com)

二、一种自创的简单、近似算法

应该有很多这样使用的,只是介绍的少,属于私下处理算法。

简单来说,思路是:

考查原始信号时间点和重采样后信号的时间点,一般来说,重采样信号的一个时间点,是处在原始信号的两个信号点之间(若重合,则直接不用计算),根据此时间点到原信号两个信号点的距离,用原信号这两个信号点插值得到重采样信号此点的信号幅值。

所以很明显,这种方法的优点是,不用考虑抽取、插值、滤波这些麻烦事,只考虑插值即可。

但是,缺点也有,就是不能统一处理,必须每个点都分别处理。

时间关系,后续,可以再测试下。

三、查到一个专利,但是算法还是有问题,未程序实现:

CN202010258902一种数字信号的任意重采样方法及系统

无法添加附件,是个问题。

主要步骤有:

简单说是用sinc函数实现。但是还有较多看不懂的。感觉不理解0052到0058的用意是啥。但是貌似,可以不管其他,直接用最后一个公式算就行了......

后续试试程序仿真...

有关sinc函数(也较辛格函数):20211003:数字滤波器前置知识,sinc函数与Sa函数_sinc函数与sa函数区别-CSDN博客

程序仿真来了,但是结果好像不对,不知道是这理论有问题还是哪里理解不对,欢迎交流了:

clc; clear all; close all;

A = 1.5;
B = 1;
f1 = 50;
f2 = 100;
Fs = 1000;
t = 0:1/Fs:1;
sig = A * cos(2 * pi *f1 * t) + B * sin(2 * pi *f2 * t); % 原始数据

M = 4; % 上采样率
N = 3;  % Fs1/Fs = M/N
INF_L = 50;

N1 = length(sig);
N2 = floor(N1*M/N);
sig_out = zeros(1,N2);

for m=1:N2
   if(m==800)
       m=800;
   end
    seq_begin = floor(M*m/N) - INF_L + 1;
    if(seq_begin<1)
        seq_begin=1;
    end
    seq_end = floor(M*m/N) + INF_L;
    if(seq_end>N1)
        seq_end=N1;
    end
    for n=seq_begin:seq_end
        tmp = abs(M*m - N*n)+1;
        sig_out(m) = sig_out(m) + sig(n)*sin(pi*tmp/M)./(pi*tmp/M);
    end    
end

sig_out2 = resample(sig,4,3);
subplot(3,1,1);
stem(sig,'MarkerFaceColor',[1 0 0]);xlim([1 100]);
title('原数据');
 
subplot(3,1,2);
stem(sig_out,'MarkerFaceColor',[1 1 0]);xlim([1 100]);
title('sinc重采样数据');

subplot(3,1,3);
stem(sig_out2,'MarkerFaceColor',[1 0 1]);xlim([1 100]);
title('resample重采样数据');


zhh = 1;


结果:

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/1411287.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

WordPress反垃圾评论插件Akismet有什么用?如何使用Akismet插件?

每次我们成功搭建好WordPress网站后&#xff0c;都可以在后台 >> 插件 >> 已安装的插件&#xff0c;在插件列表中可以看到有一个“Akismet反垃圾邮件&#xff1a;垃圾邮件保护”的插件&#xff08;个人觉得是翻译错误&#xff0c;应该是反垃圾评论&#xff09;。具…

【新书推荐】3.4 浮点型

本节必须掌握的知识点&#xff1a; 示例九 代码分析 汇编解析 浮点数的输出精度 【补充内容】 3.4.1 示例九 浮点型分为&#xff1a;单精度float、双精度double、长双精度long double。 类型 存储大小 值范围 精度 单精度 float 4字节 【1.2E-38~ 3.4E38】 6位小数 …

SpringMVC-对静态资源的访问

1.工程中加入静态资源 在webapp下创建static文件夹&#xff0c;此文件夹专门放入静态资源 2.使项目可以处理静态资源的请求 在SpringMVC配置文件中添加以下语句 1.引入命名空间 xmlns:mvc"http://www.springframework.org/schema/mvc" xsi:schemaLocation“http…

HarmonyOS鸿蒙学习笔记(23)监听Wifi状态变化

监听Wifi状态变化 前言创建接收状态变化的Bean对象创建订阅者和订阅事件参考资料&#xff1a; 前言 本篇博文通过动态订阅公共事件来说明怎么使用HarmonyOS监听Wifi状态的变化。关于动态订阅公共事件的概念&#xff0c;官网有详细说明&#xff0c;再次就不在赘述。博文相关项目…

UE5 Chaos系统 学习笔记

记得开插件&#xff1a; 1、锚点场 在锚点场范围内的物体静止且不被其他力场损坏 2、ClusterStrain 破裂效果的力 3、DisableField chaos破裂后的模拟物理在绿色范围内禁止模拟物理 4、ForceAndStrain 破裂效果的力 5、ForceAndStrainFallOff 破裂效果的力&#xff0c;但是…

代码随想录算法训练营第十一天 | 二叉树基础

代码随想录算法训练营第十一天 | 二叉树基础 文章目录 代码随想录算法训练营第十一天 | 二叉树基础1 二叉树的理论基础1.1 二叉树的类型1.2 二叉树的存储方式1.3 二叉树的遍历方式1.4 二叉树的定义 2 二叉树的递归遍历2.1 前序遍历2.2 中序遍历2.3 后序遍历 3 二叉树的迭代遍历…

【QT+QGIS跨平台编译】之十:【libbz2+Qt跨平台编译】(一套代码、一套框架,跨平台编译)

文章目录 一、libbz2介绍二、文件下载三、文件分析四、pro文件五、编译实践一、libbz2介绍 bzip2是一个基于Burrows-Wheeler 变换的无损压缩软件,压缩效果比传统的LZ77/LZ78压缩算法来得好。它是一款免费软件。可以自由分发免费使用。 bzip2能够进行高质量的数据压缩。它利用…

SpringBoot-yml文件的配置与读取

配置 值前边必须要有空格&#xff0c;作为分隔符 使用空格作为缩进表示层级关系&#xff0c;相同的层级左侧对齐 获取 使用Value(”${键名}”) 使用ConfigurationProperties(prefix "前缀") 1.前缀要与yml配置文件中的前缀一致 2.实体类的字段名与配置文件中的键名一…

linux conda 配置 stable video diffusion

安装教程 1 下载仓库源码 git clone https://github.com/Stability-AI/generative-models.git2 创建conda环境 conda create -n svd python3.10 conda activate svd3 安装pytorch gpu cuda和cudnn请参考其他链接配置&#xff0c;使用 conda 或者 pip 安装 pytorch # 使用c…

第二集《闻法仪轨》

请大家打开讲义第三面&#xff0c;甲二、于法、法师发起承事。 我们身为一个大乘的佛弟子&#xff0c;我们这一念明了的心&#xff0c;在一生当中&#xff0c;会遇到很多很多的佛法&#xff0c;也会遇到很多很多的法师&#xff0c;但不是所有的法师跟佛法对我们都是帮助的&…

EasyCVR视频智能监管系统方案设计与应用

随着科技的发展&#xff0c;视频监控平台在各个领域的应用越来越广泛。然而&#xff0c;当前的视频监控平台仍存在一些问题&#xff0c;如视频质量不高、监控范围有限、智能化程度不够等。这些问题不仅影响了监控效果&#xff0c;也制约了视频监控平台的发展。 为了解决这些问…

HCIA——29HTTP、万维网、HTML、PPP、ICMP;万维网的工作过程;HTTP 的特点HTTP 的报文结构的选择、解答

学习目标&#xff1a; 计算机网络 1.掌握计算机网络的基本概念、基本原理和基本方法。 2.掌握计算机网络的体系结构和典型网络协议&#xff0c;了解典型网络设备的组成和特点&#xff0c;理解典型网络设备的工作原理。 3.能够运用计算机网络的基本概念、基本原理和基本方法进行…

Javaweb之SpringBootWeb案例之阿里云OSS服务入门的详细解析

2.3.2 入门 阿里云oss 对象存储服务的准备工作我们已经完成了&#xff0c;接下来我们就来完成第二步操作&#xff1a;参照官方所提供的sdk示例来编写入门程序。 首先我们需要来打开阿里云OSS的官方文档&#xff0c;在官方文档中找到 SDK 的示例代码&#xff1a; 参照官方提供…

爬取第一试卷网高三数学试卷并下载到本地

import requests import re import os filename 试卷\\ if not os.path.exists(filename):os.mkdir(filename) url https://www.shijuan1.com/a/sjsxg3/list_727_1.html headers {"User-Agent":"Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.…

【QT+QGIS跨平台编译】之九:【LZ4+Qt跨平台编译】(一套代码、一套框架,跨平台编译)

文章目录 一、LZ4介绍二、文件下载三、文件分析四、pro文件五、编译实践一、LZ4介绍 LZ4是一种无损压缩算法,压缩速度为每核心400MB/s。 LZ4是目前效率最高的压缩算法,更加侧重于压缩/解压缩速度,压缩比并不突出,本质上就是时间换空间。 LZ4库是使用BSD许可证作为开放源码…

Dockerfile里ADD * 保留原来的目录结构

1、问题 给新模块写Dockerfile&#xff0c;很多静态资源分散在各个目录&#xff0c;于是Dockerfile里我直接一句&#xff1a; ADD ./* /dest/镜像出来后&#xff0c;启动容器&#xff0c;进入容器种后发现&#xff1a;文件拷贝成功&#xff0c;但原来的目录结构都不在了&…

HCIE之BGP正则表达式(四)

BGP 一、AS-Path正则表达式数字| 等同于或的关系[]和.$ 一个字符串的结束_代表任意^一个字符串的开始()括号包围的是一个组合\ 转义字符* 零个或多个&#xff1f;零个或一个一个或多个 二、BGP对等体组 一、AS-Path正则表达式 正则表达式是按照一定模版匹配字符串的公式 AR3上…

Java面试题(6)

28.创建线程池有哪几种方式 newFixedThreadPool(int nThreads) &#xff1a;创建一个固定长度的线程池&#xff0c;如果有线程发生错误而结束&#xff0c; 线程池会补充一个新线程。 newCachedThreadPool() &#xff1a;创建一个可缓存的线程池&#xff0c;会自动回收和创建空…

小黑艰难的前端啃bug之路:内联元素之间的间隙问题

今天开始学习前端项目&#xff0c;遇到了一个Bug调了好久&#xff0c;即使margin为0&#xff0c;但还是有空格。 小黑整理&#xff0c;用四种方法解决了空白问题 <!DOCTYPE html> <html><head><meta charset"utf-8"><title></tit…

关于达梦认证DCA DCP,TIDB认证PCTA PCTP考试那点事儿

文章最后有彩蛋&#xff0c;一定要看到最后... 一、正确的道路上遇到正确的你 伴随中国数据库领域的快速技术进步&#xff0c;国内数据库生态蓬勃发展&#xff0c;并不断涌现出极具创新力的产品&#xff0c;推动了数据库应用的遍地开花。截至2024年1月&#xff0c;墨天轮数据社…