C语言IIR双向滤波

设计一个0.5~1Hz的IIR滤波器，用巴特沃斯或者契比雪夫2，看零极点图是稳定的。

设计如下：

function Hd = iir_highpass_05_1_buter
%IIR_HIGHPASS_05_1_BUTER Returns a discrete-time filter object.% MATLAB Code
% Generated by MATLAB(R) 9.5 and DSP System Toolbox 9.7.
% Generated on: 27-Nov-2023 20:11:44% Butterworth Highpass filter designed using FDESIGN.HIGHPASS.% All frequency values are in Hz.
Fs = 250;  % Sampling FrequencyFstop = 0.5;         % Stopband Frequency
Fpass = 1;           % Passband Frequency
Astop = 30;          % Stopband Attenuation (dB)
Apass = 1;           % Passband Ripple (dB)
match = 'stopband';  % Band to match exactly% Construct an FDESIGN object and call its BUTTER method.
h  = fdesign.highpass(Fstop, Fpass, Astop, Apass, Fs);
Hd = design(h, 'butter', 'MatchExactly', match);
% % 
% % Get the transfer function values.
% [b, a] = tf(Hd);
% b 
% a
% close all
% figure
% freqz(b,a);
% figure
% zplane(b,a);
% [H,w]=freqz(b,a);
% % [H,w]=freqz(b,a,5000);
% N = length(H);
% Hf=abs(H);  %取幅度值实部
% Hx=angle(H);  %取相位值对应相位角
% % clf
% figure
% % plot(w,20*log10(Hf))  %幅值变换为分贝单位
% plot([0:1/N:1-1/N]*Fs/2,20*log10(Hf))  %幅值变换为分贝单位
% title('离散系统幅频特性曲线')
% figure
% % plot(w,Hx)
% plot([0:1/N:1-1/N]*Fs/2,Hx);
% title('离散系统相频特性曲线')% [EOF]

或：

function Hd = iir_highpass_05_1_cheby2
%IIR_HIGHPASS_05_1_BUTER Returns a discrete-time filter object.% MATLAB Code
% Generated by MATLAB(R) 9.5 and DSP System Toolbox 9.7.
% Generated on: 27-Nov-2023 20:10:02% Chebyshev Type II Highpass filter designed using FDESIGN.HIGHPASS.% All frequency values are in Hz.
Fs = 250;  % Sampling FrequencyFstop = 0.5;         % Stopband Frequency
Fpass = 1;           % Passband Frequency
Astop = 30;          % Stopband Attenuation (dB)
Apass = 1;           % Passband Ripple (dB)
match = 'stopband';  % Band to match exactly% Construct an FDESIGN object and call its CHEBY2 method.
h  = fdesign.highpass(Fstop, Fpass, Astop, Apass, Fs);
Hd = design(h, 'cheby2', 'MatchExactly', match);
% 
% % Get the transfer function values.
% [b, a] = tf(Hd);
% b 
% a
% close all
% figure
% freqz(b,a);
% figure
% zplane(b,a);
% % [H,w]=freqz(b,a);
% [H,w]=freqz(b,a,5000);
% N = length(H);
% Hf=abs(H);  %取幅度值实部
% Hx=angle(H);  %取相位值对应相位角
% % clf
% figure
% % plot(w,20*log10(Hf))  %幅值变换为分贝单位
% plot([0:1/N:1-1/N]*Fs/2,20*log10(Hf))  %幅值变换为分贝单位
% title('离散系统幅频特性曲线')
% figure
% % plot(w,Hx)
% plot([0:1/N:1-1/N]*Fs/2,Hx);
% title('离散系统相频特性曲线')% [EOF]

可以得到滤波器系数a,b.

设计一个通用的双向滤波程序，包括滤波函数和数组逆序函数,包含可以和matlab对数据的写数据代码。主要有：

main.c

#include"stdio.h"
#include"stdlib.h"
#include"math.h"#include"highpass_iir_05_1.h"#define M_PI 3.1415
#define N 200
#define F1 3
#define F2 40//float fir_coef_b[5] = { 0.979965406517890, - 3.91970687664644,	5.87948294331178, - 3.91970687664644,	0.979965406517890 };
//float fir_coef_a[5] = { 1, - 3.95936787226356,	5.87908160514721, - 3.88004583425791,	0.960332197971774 };float fir_coef_b[5] = { 1,2,1,3,1 };
float fir_coef_a[5] = { 1,1,1,1,1 };int  main()
{// 生成正弦信号float a[N] = { 0 }, b[N] = { 0 };for (int i = 0; i < N; i++){a[i] = (float)(sin(2 * M_PI * i * F1 / N) + 0.25 * sin(2 * M_PI * i * F2 / N));}FILE* fp;fp = fopen("text.txt", "w");if (feof(fp)){printf("NULL");exit(0);//表示如果读取为空文件就正常退出}for (int i = 0; i < N; i++)fprintf(fp, "%.6f\n", a[i]);fclose(fp);int m;// 双向滤波m = iir_filter_zhh(a, b, N, fir_coef_a, 5, fir_coef_b,5);reverse_arry_float(b, N);memcpy(a, b, N * sizeof(float));m = iir_filter_zhh(a, b, N, fir_coef_a, 5, fir_coef_b, 5);reverse_arry_float(b, N);fp = fopen("text_fir.txt", "w");if (feof(fp)){printf("NULL");exit(0);//表示如果读取为空文件就正常退出}for (int i = 0; i < N; i++)fprintf(fp, "%.6f\n", b[i]);fclose(fp);return 0;
}

highpass_iir_05_1.c

#include"stdio.h"
#include"stdlib.h"
#include"string.h"#include "assert.h" // MATLAB函数
//function y = iir_filter_zhh(b, a, x)
//
//len_x = length(x);
//y = zeros(size(x));
//
//len_b = length(b);
//len_a = length(a);
//x_buf = [zeros(len_b - 1, 1); x];
//y_buf = [zeros(len_a - 1, 1); y];
//
//b_fir = fliplr(b);
//a_fir = fliplr(a);
//
//for i = 1:len_x
//y(i) = 1 / a_fir(end) * (b_fir * x_buf(i:i + len_b - 1) - a_fir(1:end - 1) * y_buf(i:i + len_a - 2));
//y_buf(i + len_a - 1) = y(i);
//end
//
//endint iir_filter_zhh(float* sig_in, float* sig_out, int sig_len, float* a, int a_len, float* b, int b_len)
{// 参数个数检查,至少3个至多5个参数.系统自行检查，参数个数不对及类型不对会自动报错。assert(sig_len > 1);assert(a_len > 1);assert(b_len > 1);int i, j;float tmp1, tmp2;float* sig_in_buf = (float*)malloc((size_t)(sig_len + b_len - 1) * sizeof(float));if (sig_in_buf == NULL) {//判空printf("sig_in_buf malloc error");//打印错误信息return 1;}memset(sig_in_buf, 0, (size_t)(b_len - 1) * sizeof(float));memcpy(sig_in_buf + b_len - 1, sig_in, (size_t)sig_len * sizeof(float));float* sig_out_buf = (float*)malloc((size_t)(sig_len + a_len - 1) * sizeof(float));if (sig_out_buf == NULL) {//判空printf("sig_out_buf malloc error");//打印错误信息return 1;}memset(sig_out_buf, 0, (size_t)(sig_len + a_len - 1) * sizeof(float));for (i = 0; i < sig_len; i++) {tmp1 = 0; tmp2 = 0;for (j = 0; j < b_len; j++) {tmp1 += sig_in_buf[i - j + b_len - 1] * b[j];}for (j = 1; j < a_len; j++) {tmp2 += sig_out_buf[i - j + a_len - 1] * a[j];}sig_out[i] = (tmp1 - tmp2) / a[0];sig_out_buf[i + a_len - 1] = sig_out[i];}free(sig_in_buf);free(sig_out_buf);return 0;
}// 数组逆序 float版
void reverse_arry_float(float* arry, unsigned arry_len) // arry_len可以为1，不报错，此时arry只有一个元素，调用后不变
{int i = 0;  //循环变量1, i的值为数组第一个元素的下标int j = arry_len - 1;  //循环变量2, j的值为数组最后一个元素的下标float fir_idx_buf;  //互换时的中间存储变量for (; i < j; ++i, --j)  /*因为i和j已经初始化过了, 所以表达式1可以省略, 但表达式1后面的分号不能省。*/{fir_idx_buf = arry[i];arry[i] = arry[j];arry[j] = fir_idx_buf;}
}

highpass_iir_05_1.h

#pragma onceextern int iir_filter_zhh(float* sig_in, float* sig_out, int sig_len, float* a, int a_len, float* b, int b_len);
extern void reverse_arry_float(float* arry, unsigned arry_len);

matlab对数据代码：

close all,clear,clc% fs = 200;
% f1 = 3;f2 = 40;
% t = 0:1/fs:1-1/fs;
% x = sin(2*pi*t*f1)+0.25*sin(2*pi*t*f2);fid=fopen('D:\zhh\work\VS projects\highpass_iir_05_1\text.txt'); %D:\zhh\work\VS projects\filter_int_v3
x=fscanf(fid,'%f');
fclose(fid);
filename = 'data.txt';
% dlmwrite(filename, x);
fid = fopen(filename, 'w');
for i=1:length(x)fprintf(fid, '%.6f\n', x(i));
end
fclose(fid);% Hd = iir_highpass_05_1_cheby2;
% % 直接matlab滤波
% [b, a] = tf(Hd);
b=[1,2,1,3,1 ];
a=[1,1,1,1,1 ];
y1 = iir_filter_zhh(b,a,x);
y2 = flipud(  filter( b,a,flipud(y1) )  );figure
plot(x);hold on
plot(y1);
figure
plot(x);hold on
plot(y2);filename = 'data_fir.txt';
fid = fopen(filename, 'w');
for i=1:length(x)fprintf(fid, '%.6f\n', y2(i));
end
fclose(fid);zhh = 1;

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

如若内容造成侵权/违法违规/事实不符，请联系多彩编程网进行投诉反馈email:809451989@qq.com，一经查实，立即删除！

相关文章

流畅的Python （节选）

流畅的Python （节选）

0 前言 节选学习部分有用的内容 Fluent Python 2 序列构成的数组 Python 会忽略代码里 []、{} 和 () 中的换行&#xff0c;因此如果你的代码里有多行的列表、列表推导、生成器表达式、字典这一类的&#xff0c;可以省略不太好看的续行符 \。 元组其实是对数据的记录&#x…
阅读更多...
苹果输入法怎么换行？3个换行技巧，速速掌握！

苹果输入法怎么换行？3个换行技巧，速速掌握！

在日常打字的时候&#xff0c;我们经常会遇到需要换行的情况。比如&#xff0c;在聊天、写作、编辑文档等场景下&#xff0c;当一行文字输入完成后&#xff0c;我们通常需要将光标移动到下一行再继续输入文字。那么这时候就需要我们进行换行操作。 然而&#xff0c;很多用户对…
阅读更多...
webpack开发环境

webpack开发环境

文章目录 前言webpack.config.js使用 source mapwebpack.config.jssrc/print.js 选择一个开发工具使用 watch mode(观察模式)package.jsonsrc/print.js使用 webpack-dev-serverwebpack.config.js package.json使用 webpack-dev-middlewareprojectserver.jspackage.jsonpackage.…
阅读更多...
LeetCode | 二叉树的最大深度

LeetCode | 二叉树的最大深度

LeetCode | 二叉树的最大深度 OJ链接 这里需要注意的一点是每次有返回值&#xff0c;需要定义变量来保存上一次的值最后取最高的一方加1 int maxDepth(struct TreeNode* root) {if(root NULL)return NULL;int left maxDepth(root->left);int right maxDepth(root->r…
阅读更多...
秒杀业务

秒杀业务

1.缓存秒杀商品库存到redis 2.用户访问秒杀商品&#xff0c;进行商品的抢购 3.判断用户是否存在未支付的秒杀订单&#xff0c;如果存在&#xff0c;告知用户&#xff0c;请先支付。&#xff08;幂等&#xff09; 4.判断排队人数是否达到上线。redis ,incrment() 达到上线&am…
阅读更多...
【广州华锐视点】3D宪法普法知识宣传展厅——线上法律知识学习新途径

【广州华锐视点】3D宪法普法知识宣传展厅——线上法律知识学习新途径

随着科技的不断发展&#xff0c;人们的生活方式也在不断地改变。在这个信息爆炸的时代&#xff0c;传统的普法教育方式已经无法满足人们的需求。为了适应这一变化&#xff0c;越来越多的教育机构开始尝试利用现代科技手段进行普法教育。其中&#xff0c;3D宪法普法知识宣传展厅…
阅读更多...
HarmonyOS——UI开展前的阶段总结

HarmonyOS——UI开展前的阶段总结

当足够的了解了HarmonyOS的相关特性之后&#xff0c;再去介入UI&#xff0c;你会发现无比的轻松&#xff0c;特别当你有着其他的声明式UI开发的经验时&#xff0c;对于HarmonyOS的UI&#xff0c;大致一扫&#xff0c;也就会了。 如何把UI阐述的简单易懂&#xff0c;又能方便大…
阅读更多...
【Java】3. 字面量

【Java】3. 字面量

3.字面量 字面量类型说明程序中的写法整数不带小数的数字666&#xff0c;-88小数带小数的数字13.14&#xff0c;-5.21字符必须使用单引号&#xff0c;有且仅能一个字符‘A’&#xff0c;‘0’&#xff0c; ‘我’字符串必须使用双引号&#xff0c;内容可有可无“HelloWorld”&…
阅读更多...
[论文阅读]CT3D——逐通道transformer改进3D目标检测

[论文阅读]CT3D——逐通道transformer改进3D目标检测

CT3D 论文网址&#xff1a;CT3D 论文代码&#xff1a;CT3D 简读论文 本篇论文提出了一个新的两阶段3D目标检测框架CT3D,主要的创新点和方法总结如下: 创新点: (1) 提出了一种通道注意力解码模块,可以进行全局和局部通道聚合,生成更有效的解码权重。 (2) 提出了建议到点嵌…
阅读更多...
235. 二叉搜索树的最近公共祖先

235. 二叉搜索树的最近公共祖先

235. 二叉搜索树的最近公共祖先 原题链接&#xff1a;完成情况&#xff1a;解题思路&#xff1a;_235二叉搜索树的最近公共祖先_递归_235二叉搜索树的最近公共祖先_简洁递归 参考代码&#xff1a;错误经验吸取 原题链接&#xff1a; 235. 二叉搜索树的最近公共祖先 https://…
阅读更多...
zookeeper集群+kafka集群

zookeeper集群+kafka集群

Kafka3.0之前依赖于zookeeper Zookeeper开源&#xff0c;分布式的架构&#xff0c;提供协调服务&#xff08;apache项目&#xff09; 基于观察者模式设计的分布式服务管理架构 存储和管理数据&#xff0c;分布式节点上的服务接受观察者的注册&#xff0c;一旦分布式节点上的…
阅读更多...
ubuntu22安装vmtool正确姿势

ubuntu22安装vmtool正确姿势

确保GUI处于X11 sudo vi /etc/gdm3/custom.conf取消#WaylandEnablefalse的注释systemctl restart gdm3#重启gdm3屏幕会空白一分钟&#xff0c;但不要担心&#xff0c;当它恢复时&#xff0c;GUI应该在X11模式 安装openvmtools 确保包索引已更新: sudo apt-get updatesudo ap…
阅读更多...
AI聊天 AI绘画 AI视频 AI制作PPT

AI聊天 AI绘画 AI视频 AI制作PPT

文章目录&#xff1a; 一&#xff1a;AI聊天 二&#xff1a;AI绘画 三&#xff1a;AI视频 四&#xff1a;AI制作PPT 这里主要放一些国内我感觉好用的&#xff0c;国外或者更多请移步——>AI-Chat_Draw_Video_PPT 一&#xff1a;AI聊天 文心一言&#xff1a;百度旗下&a…
阅读更多...
在Linux本地部署开源自托管导航页配置服务Dashy并远程访问

在Linux本地部署开源自托管导航页配置服务Dashy并远程访问

文章目录 简介1. 安装Dashy2. 安装cpolar3.配置公网访问地址4. 固定域名访问 简介 Dashy 是一个开源的自托管的导航页配置服务&#xff0c;具有易于使用的可视化编辑器、状态检查、小工具和主题等功能。你可以将自己常用的一些网站聚合起来放在一起&#xff0c;形成自己的导航…
阅读更多...
Spring 日志

Spring 日志

日志的作用: 1.定位和发现问题 2.系统监控 3.数据采集 观察日志 先写一段打印日志的代码 日志内容 日志级别分类 默认日志级别是Info,级别一下的就不打印了 Spring 帮我们集成了日志框架,我们直接使用即可 我们测试一下用日志框架打印日志是如何 我们就会发现打印的结果跟…
阅读更多...
15 动态规划解统计全为1的正方形子矩阵

15 动态规划解统计全为1的正方形子矩阵

来源&#xff1a;LeetCode第1277题 难度&#xff1a;中等 描述&#xff0c;给你一个m*n的矩阵&#xff0c;矩阵中的元素不是0就是1&#xff0c;请你统计并返回其中完全由1组成的正方形子矩形的个数&#xff1b; 分析&#xff1a;可以使用动态规划求解dp[i][j]表示以[i][j]为…
阅读更多...
vue实现el-tree操作后默认展开该节点和同级节点拖拽并进行前后端排序

vue实现el-tree操作后默认展开该节点和同级节点拖拽并进行前后端排序

问题一&#xff1a;实现el-tree 删除、添加、编辑后默认展开该节点 操作视频如下 el-tree节点操作后仍展开 节点代码 <template><el-treev-loading"loading"ref"tree"element-loading-text"加载中"highlight-current:data"treeD…
阅读更多...
【Web系列二十七】Vue实现dom元素拖拽并限制移动范围

【Web系列二十七】Vue实现dom元素拖拽并限制移动范围

目录 需求 拖拽功能封装 使用拖拽功能 vite-env.d.ts main.ts test.vue 需求 dom元素拖拽并限制在父组件范围内 拖拽功能封装 export const initVDrag (vue) > {vue.directive(drag, (el) > {const oDiv el // 当前元素oDiv.onmousedown (e) > {let target…
阅读更多...
Spark---创建DataFrame的方式

Spark---创建DataFrame的方式

1、读取json格式的文件创建DataFrame 注意&#xff1a; 1、可以两种方式读取json格式的文件。 2、df.show()默认显示前20行数据。 3、DataFrame原生API可以操作DataFrame。 4、注册成临时表时&#xff0c;表中的列默认按ascii顺序显示列。 df.createTempView("mytab…
阅读更多...
《HelloGitHub》第 92 期

《HelloGitHub》第 92 期

兴趣是最好的老师&#xff0c;HelloGitHub 让你对编程感兴趣&#xff01; 简介 HelloGitHub 分享 GitHub 上有趣、入门级的开源项目。 https://github.com/521xueweihan/HelloGitHub 这里有实战项目、入门教程、黑科技、开源书籍、大厂开源项目等&#xff0c;涵盖多种编程语言 …
阅读更多...
最新文章

Copyright @ 2022~2023