① 二阶滤波器用c语言怎么写
这个可比你想象的复杂多了,s是个复变量,1/(s+1)极点在-1,要想用C语言写,必须理解清楚下面几个问题:
1、输入必须是个有限序列,比如(x+yi),x和y分别是两个长度为N的数组
2、要过滤的频率,必须是个整型值,或者是个整型区间
3、输出结果同样是两个长度为N的数组(p+qi)
4、整个程序需要使用最基本的复数运算,这一点C语言本身不提供,必须手工写复函数运算库
5、实现的时候具体算法还需要编,这里才是你问题的核心。
我可以送你一段FFT的程序,自己琢磨吧,和MATLAB的概念差别很大:
#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <windows.h>
#include "complex.h"
extern "C" {
// Discrete Fourier Transform (Basic Version, Without Any Enhancement)
// return - Without Special Meaning, constantly, zero
int DFT (long count, CComplex * input, CComplex * output)
{
assert(count);
assert(input);
assert(output);
CComplex F, X, T, W; int n, i;
long N = abs(count); long Inversing = count < 0? 1: -1;
for(n = 0; n < N ; n++){ // compute from line 0 to N-1
F = CComplex(0.0f, 0.0f); // clear a line
for(i = 0; i < N; i++) {
T = input[i];
W = HarmonicPI2(Inversing * n * i, N);
X = T * W;
F += X; // fininshing a line
}//next i
// save data to outpus
memcpy(output + n, &F, sizeof(F));
}//next n
return 0;
}//end DFT
int fft (long count, CComplex * input, CComplex * output)
{
assert(count);
assert(input);
assert(output);
int N = abs(count); long Inversing = count < 0? -1: 1;
if (N % 2 || N < 5) return DFT(count, input, output);
long N2 = N / 2;
CComplex * iEven = new CComplex[N2]; memset(iEven, 0, sizeof(CComplex) * N2);
CComplex * oEven = new CComplex[N2]; memset(oEven, 0, sizeof(CComplex) * N2);
CComplex * iOdd = new CComplex[N2]; memset(iOdd , 0, sizeof(CComplex) * N2);
CComplex * oOdd = new CComplex[N2]; memset(oOdd , 0, sizeof(CComplex) * N2);
int i = 0; CComplex W;
for(i = 0; i < N2; i++) {
iEven[i] = input[i * 2];
iOdd [i] = input[i * 2 + 1];
}//next i
fft(N2 * Inversing, iEven, oEven);
fft(N2 * Inversing, iOdd, oOdd );
for(i = 0; i < N2; i++) {
W = HarmonicPI2(Inversing * (- i), N);
output[i] = oEven[i] + W * oOdd[i];
output[i + N2] = oEven[i] - W * oOdd[i];
}//next i
return 0;
}//end FFT
void __stdcall FFT(
long N, // Serial Length, N > 0 for DFT, N < 0 for iDFT - inversed Discrete Fourier Transform
double * inputReal, double * inputImaginary, // inputs
double * AmplitudeFrequences, double * PhaseFrequences) // outputs
{
if (N == 0) return;
if (!inputReal && !inputImaginary) return;
short n = abs(N);
CComplex * input = new CComplex[n]; memset(input, 0, sizeof(CComplex) * n);
CComplex * output= new CComplex[n]; memset(output,0, sizeof(CComplex) * n);
double rl = 0.0f, im = 0.0f; int i = 0;
for (i = 0; i < n; i++) {
rl = 0.0f; im = 0.0f;
if (inputReal) rl = inputReal[i];
if (inputImaginary) im = inputImaginary[i];
input[i] = CComplex(rl, im);
}//next i
int f = fft(N, input, output);
double factor = n;
//factor = sqrt(factor);
if (N > 0)
factor = 1.0f;
else
factor = 1.0f / factor;
//end if
for (i = 0; i < n; i++) {
if (AmplitudeFrequences) AmplitudeFrequences[i] = output[i].getReal() * factor;
if (PhaseFrequences) PhaseFrequences[i] = output[i].getImaginary() * factor;
}//next i
delete [] output;
delete [] input;
return ;
}//end FFT
int __cdecl main(int argc, char * argv[])
{
fprintf(stderr, "%s usage:\n", argv[0]);
fprintf(stderr, "Public Declare Sub FFT Lib \"wfft.exe\" \
(ByVal N As Long, ByRef inputReal As Double, ByRef inputImaginary As Double, \
ByRef freqAmplitude As Double, ByRef freqPhase As Double)");
return 0;
}//end main
};//end extern "C"
② 巴特沃斯滤波器c语言实现
说的很对,滤波玩的就是增益(衰减)变化,不同的频率,不同的增益(衰减)。称幅频曲线。
1、巴特奥斯滤波器的截止频率指-3db通频带频率,也就是在这个频率以内保证畅通(通带)。
2、另一个指标叫做阻带,频率大于此值能够保证衰减大于某值
③ 求IIR及FIR数字滤波器的C语言实现。(VC++)
这个问题比较复杂,最近本人也在研究数字滤波,
结合图片说一下
第一个图是fir的流程图,其中Z-1是延迟,是单个采样时间1/fs
n阶的fir滤波器就是选取最近的n+1个样本,然后使他们各自乘以自己的滤波器系数即图中的F(n),[一般其他书的表示是h(n)]
然后相加得到输出的y(n)就是一个输出点
,其中F(n)的得出需要根据采样频率和滤波器的通带和阻带来决定
其中为了改善旁瓣的幅值,一般在采样后给样本或者h(n)加窗,当然可以用“最佳方法”来做
得出h(n)大致方法是先将矩形窗进行DFT,得出h(n),然后对h(n)进行加窗得出h(k),然后将∑h(k)×x(n)=y(n),假如阶数较多可以用傅里叶变换使时域变频域后再将卷积相加,可以利用FFT来改进实时性,提升速度
上面就是fir滤波器的简述
第二个图片上传不了,直接给链接
http://image..com/i?ct=503316480&z=0&tn=imagedetail&word=%D2%BB%BD%D7iir%C2%CB%B2%A8%C6%F7&in=12708&cl=2&cm=1&sc=0&lm=-1&pn=0&rn=1&di=2607528304&ln=1054&fr=
图中的Z-1是延时,iir滤波器也叫无限冲击响应滤波器,是有反馈的,
图中的是一阶的,相对fir滤波器来说,iir滤波器可以用较低的阶数来获得较好的滤波特效。但是其相位特性较差。
鉴于实用性,还是建议楼主去图书馆借书看,网络不可能得到确实的方案,
楼主可以去借“数字信号处理”的书,国外的中译本就有详细介绍fir和iir以及fft还有其他变换,国内的dsp大都几乎是dsp用户手册的中译本,对上述问题都是很简陋地带过,不予置评。
本人推荐一本书在www.ouravr.com上面的dsp专栏有下载,40多M,叫DSP算法、应用和设计,本人有这本实体书,写的较好
④ c语言中值滤波问题
1. 是规定做中值滤波的点不含边缘的点(取决于中值滤波窗口大小)。 2,对图像边缘部分的信息进行镜像处理。
⑤ 帮帮忙,能不能给我 基于C语言的FIR滤波器设计的程序代码(包括CMD,C,ASM),谢谢了 真的很急!!!
#include"math.h"
void firwin(n,band,fln,fhn,wn,h)
int n,band,wn;
double fln,fhn,h[];
{int i,n2,mid;
double s,pi,wc1,wc2,beta,delay;
double window();
beta=0.0;
if(wn==7)
{printf("input beta parameter of Kaiser window(2<beta<10)\n");
scanf("%1f",&beta);
}
pi=4.0*atan(1.0);
if((n%2)==0)/*如果n是偶数*/
{n2=n/2+1;/*这行什么意思*/
mid=1;
}
else
{n2=n/2;
mid=0;
}
delay=n/2.0;
wc1=2.0*pi*fln;
if(band>=3) wc2=2.0*pi*fhn;/*先判断用户输入的数据,如果band参数大于3*/
switch(band)
{case 1:
{for(i=0;i<=n2;i++)
{s=i-delay;
h[i]=(sin(wc1*s)/(pi*s))*window(wn,n+1,i,beta);
h[n-i]=h[i];
}
if(mid==1) h[n/2]=wc1/pi;
break;
}
case 2:
{for(i=0;i<=n2;i++)
{s=i-delay;
h[i]=(sin(pi*s)-sin(wc1*s))/(pi*s);
h[i]=h[i]*window(wn,n+1,i,beta);
h[n-i]=h[i];
}
if(mid==1) h[n/2]=1.0-wc1/pi;
break;
}
case 3:
{for(i=0;i<n2;i++)
{s=i-delay;
h[i]=(sin(wc2*s)-sin(wc1*s))/(pi*s);
h[i]=h[i]*window(wn,n+1,i,beta);
h[n-i]=h[i];
}
if(mid==1)h[n/2]=(wc2-wc1)/pi;
break;
}
case 4:
{for(i=0;i<=n2;i++)
{s=i-delay;
h[i]=(sin(wc1*s)+sin(pi*s)-sin(wc2*s))/(pi*s);
h[i]=h[i]*window(wn,n+1,i,beta);
h[n-i]=h[i];
}
if(mid==1)h[n/2]=(wc1+pi-wc2)/pi;
break;
}
}
}
static double window(type,n,i,beta)
int i,n,type;
double beta;
{int k;
double pi,w;
double kaiser();
pi=4.0*atan(1.0);
w=1.0;
switch(type)
{case 1:
{w=1.0;
break;
}
case 2:
{k=(n-2)/10;
if(i<=k)
w=0.5*(1.0-cos(i*pi/(k+1)));
break;
}
case 3:
{w=1.0-fabs(1.0-2*i/(n-1.0));
break;
}
case 4:
{w=0.5*(1.0-cos(2*i*pi/(n-1)));
break;
}
case 5:
{w=0.54-0.46*cos(2*i*pi/(n-1));
break;
}
case 6:
{w=0.42-0.5*cos(2*i*pi/(n-1))+0.08*cos(4*i*pi/(n-1));
break;
}
case 7:
{w=kaiser(i,n,beta);
break;
}
}
return(w);
}
static double kaiser(i,n,beta)
int i,n;
double beta;
{
double a,w,a2,b1,b2,beta1;
double bessel0();
b1=bessel0(beta);
a=2.0*i/(double)(n-1)-1.0;
a2=a*a;
beta1=beta*sqrt(1.0-a2);
b2=bessel0(beta1);
w=b2/b1;
return(w);
}
static double bessel0(x)
double x;
{int i;
double d,y,d2,sum;
y=x/2.0;
d=1.0;
sum=1.0;
for(i=1;i<=25;i++)
{d=d*y/i;
d2=d*d;
sum=sum+d2;
if(d2<sum*(1.0e-8)) break;
}
return(sum);
}
这是窗函数法的,当然还有其他的比如切比雪夫,零相位滤波什么的,我也在研究,不是很懂哈
⑥ 给定一个信号s,包括高频、低频成份,现在做一个低通滤波器对信号进行滤波,求c语言程序,急急急!
教你实现方法吧:
第一步:你需要采用MATLAB 设计滤波器的脉冲响应序列,如记为h(n) = [ h1,h2,... hN]
第二步:c 语言实现的是 y(n) = h(n) 卷积x(n) 其中h(n)和x(n)都是向量
详细代码很简单,自己来吧。
⑦ 只是用简单的C语言,能实现高斯滤波吗
可以的。
int*** SmoothImage(int ***XImage ,int width, int height, int channel)
{
double sigma = 1.85; //(n/2 -1)*0.3 +0.8 { n = 9 ,no. of elements}
double conv[3][3];
double hg = 0;
//Convolution kernel
for(i=0; i<3;i++)
{
for(j=0;j<3;j++)
{
int u=i-1; //subtract from centre element index in 3*3 (1,1)
int v=j-1;
conv[i][j] =exp ( - ((u*u) + (v*v)) / (2*sigma*sigma) );
hg += conv[i][j];
}
}
for(int i=0; i<3;i++)
{
for(int j=0;j<3;j++)
{
conv[i][j] = conv[i][j] / hg;
}
}
int*** sXImage = 0;
sXImage = CreateImageMatrix( sXImage , width , height , channel );// allocating a 3d array
//Assigning weights
for(int i =0; i < height ; i++)
{
for(int j =0; j < width ; j++)
{
for(int k =0; k < 3 ; k++)
{
double val = 0;
double valw =0;
if(j-1 > 0 && i-1 > 0)
{
val += conv[0][0] * XImage[i-1][j-1][k];
valw += conv[0][0];
}
if(i-1 > 0 )
{
val += conv[0][1] * XImage[i-1][j][k];
valw += conv[0][1];
}
if(i-1 > 0 && j+1 < width)
{
val += conv[0][2] * XImage[i-1][j+1][k];
valw += conv[0][2];
}
if(j-1 > 0 )
{
val += conv[1][0] * XImage[i][j-1][k];
valw += conv[1][0];
}
val += conv[1][1] * XImage[i][j][k];
valw += conv[1][1];
if(j+1 < width)
{
val += conv[1][2] * XImage[i][j+1][k];
valw += conv[1][2];
}
if(j-1 > 0 && i+1 < height)
{
val += conv[2][0] * XImage[i+1][j-1][k];
valw += conv[2][0];
}
if(i+1 < height)
{
val += conv[2][1] * XImage[i+1][j][k];
valw += conv[2][1];
}
if(j+1 < width && i+1 < height)
{
val += conv[2][2] * XImage[i+1][j+1][k];
valw += conv[2][2];
}
sXImage[i][j][k] = val / valw;
}
}
}
return( sXImage);
}
⑧ 求:一个关于FIR带通滤波器的C语言设计程序 代码
short h[], short y[])
{
int i, j, sum; for (j = 0; j < 100; j++) {
sum = 0;
for (i = 0; i < 32; i++)
sum += x[i+j] * h[i];
y[j] = sum >> 15;
}
}
2
void fir(short x[], short h[], short y[])
{
int i, j, sum0, sum1;
short x0,x1,h0,h1; for (j = 0; j < 100; j+=2) {
sum0 = 0;
sum1 = 0;
x0 = x[j];
for (i = 0; i < 32; i+=2){
x1 = x[j+i+1];
h0 = h[i];
sum0 += x0 * h0;
sum1 += x1 * h0;
x0 = x[j+i+2];
h1 = h[i+1];
sum0 += x1 * h1;
sum1 += x0 * h1;
}
y[j] = sum0 >> 15;
y[j+1] = sum1 >> 15;
}
}
3
void fir(short x[], short h[], short y[])
{
int i, j, sum0, sum1;
short x0,x1,x2,x3,x4,x5,x6,x7,h0,h1,h2,h3,h4,h5,h6,h7; for (j = 0; j < 100; j+=2) {
sum0 = 0;
sum1 = 0;
x0 = x[j];
for (i = 0; i < 32; i+=8){
x1 = x[j+i+1];
h0 = h[i];
sum0 += x0 * h0;
sum1 += x1 * h0;
x2 = x[j+i+2];
h1 = h[i+1];
sum0 += x1 * h1;
sum1 += x2 * h1;
x3 = x[j+i+3];
h2 = h[i+2];
sum0 += x2 * h2;
sum1 += x3 * h2;
x4 = x[j+i+4];
h3 = h[i+3];
sum0 += x3 * h3;
sum1 += x4 * h3;
x5 = x[j+i+5];
h4 = h[i+4];
sum0 += x4 * h4;
sum1 += x5 * h4;
x6 = x[j+i+6];
h5 = h[i+5];
sum0 += x5 * h5;
sum1 += x6 * h5;
x7 = x[j+i+7];
h6 = h[i+6];
sum0 += x6 * h6;
sum1 += x7 * h6;
x0 = x[j+i+8];
h7 = h[i+7];
sum0 += x7 * h7;
sum1 += x0 * h7;
}
y[j] = sum0 >> 15;
y[j+1] = sum1 >> 15;
}
}