二阶滤波器用C语言怎么写

二阶滤波器用C语言怎么写,第1张

这个可比你想象晌大的复杂多了,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, Wint n, i

long N = abs(count)long Inversing = count <0? 1: -1

for(n = 0n <N n++){ // compute from line 0 to N-1

F = CComplex(0.0f, 0.0f)// clear a line

for(i = 0i <Ni++) {

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 = 0CComplex W

for(i = 0i <N2i++) {

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 = 0i <N2i++) {

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.0fint i = 0

for (i = 0i <ni++) {

rl = 0.0fim = 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 = 0i <ni++) {

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"

可以告诉你方法:算数平均滤波,就是求出k次采样值的总和,再除以k;中值滤波法,是把k个采样值按照从小到大排列顺序,然后找到位于最中间的那个值;最后一种差薯正不知道你们老师的防脉冲是什么意思,猜测可能是去掉k次采样值中大于虚悔或小于某个值,剩余值求平均数。手迅

让别人免费给你写程序基本上不可能,这个得花时间和精力。

可以察皮的。

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=0i<3i++)

{

for(j=0j<3j++)

{

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=0i<3i++)

{

for(int j=0j<3j++)

{

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 =0i <height i++)

{

for(int j =0j <width j++)

{

for(int k =0k <谈铅 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)

}


欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/yw/12516787.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2023-05-26
下一篇 2023-05-26

发表评论

登录后才能评论

评论列表(0条)

保存