二阶滤波器用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"

float middle_filter(float middle_value [] , intcount)

{

    float sample_value, data

    int i, j

    for (i=1 i for(j=count-1 j>=i,--j){

        if(middle_value[j-1]=middle_value[j]{

            data=middle_value[j-1]

            middle_value[j-1]=middle_value[j]

            middle_value[j]=data

        }

    }

    sample_value=middle_value(count-1)/2]

    return(sample_value)

}


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

原文地址:https://54852.com/yw/11751297.html

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

发表评论

登录后才能评论

评论列表(0条)

    保存