C语言计算数字滤波器的幅频响应和相频响应

ж•°еӯ—дҝЎеҸ·еӨ„зҗҶCиҜӯиЁҖзЁӢеәҸйӣҶпјҢP168

ж–№жі•з®Җд»Ӣпјҡ

еӯҗеҮҪж•°иҜӯеҸҘпјҡ

void gain(double b[], double a[], int m, int n, double x[], double y[], int len, int sign)

зі»ж•°д»Ӣз»Қпјҡ

bпјҡй•ҝеәҰдёәm+1пјҢеӯҳж”ҫж»ӨжіўеҷЁеҲҶеӯҗеӨҡйЎ№ејҸзҡ„зі»ж•°bпјҲiпјү
aпјҡй•ҝеәҰдёәn+1пјҢеӯҳж”ҫж»ӨжіўеҷЁеҲҶжҜҚеӨҡйЎ№ејҸзҡ„зі»ж•°aпјҲi)
m: ж»ӨжіўеҷЁеҲҶеӯҗеӨҡйЎ№ејҸзҡ„йҳ¶ж•°
n: ж»ӨжіўеҷЁеҲҶжҜҚеӨҡйЎ№ејҸзҡ„йҳ¶ж•°
x: й•ҝеәҰдёәlen,signдёә0ж—¶пјҢеӯҳж”ҫж»ӨжіўеҷЁйў‘зҺҮе“Қеә”зҡ„е®һйғЁRe[H(w)],еҪ“sign=1ж—¶пјҢеӯҳж”ҫж»ӨжіўеҷЁе№…еәҰе“Қеә”|HпјҲwпјү|пјӣеҪ“sign=2ж—¶пјҢеӯҳж”ҫз”ЁеҲҶиҙқиЎЁзӨәзҡ„ж»ӨжіўеҷЁе№…йў‘е“Қеә”|HпјҲwпјү|
yпјҡй•ҝеәҰдёәlenпјҢеҪ“sign=0ж—¶пјҢеӯҳж”ҫж»ӨжіўеҷЁйў‘зҺҮе“Қеә”зҡ„иҷҡйғЁIm[H(w)],еҪ“sign=1е’Ң2ж—¶пјҢеӯҳж”ҫж»ӨжіўеҷЁзҡ„зӣёйў‘е“Қеә”
lenпјҡйў‘зҺҮе“Қеә”зҡ„й•ҝеәҰ
signпјҡеҪ“sign=0ж—¶пјҢи®Ўз®—ж»ӨжіўеҷЁйў‘зҺҮе“Қеә”зҡ„е®һйғЁRe[H(w)]е’ҢиҷҡйғЁIm[HпјҲw)],еҪ“sign=1ж—¶пјҢи®Ўз®—ж»ӨжіўеҷЁзҡ„е№…йў‘е“Қеә”|H(w)|е’Ңзӣёйў‘е“Қеә”пјӣеҪ“sign=2ж—¶пјҢи®Ўз®—ж»ӨжіўеҷЁзҡ„е№…йў‘е“Қеә”|HпјҲwпјү|е’Ңзӣёйў‘е“Қеә”гҖӮ

еӯҗзЁӢеәҸд»Јз ҒеҰӮдёӢпјҡ

void gain(double b[], double a[], int m, int n, double x[], double y[], int len, int sign)
{
    int i, k;
    double ar, ai, br, bi, zr, zi, im, re, den, numr, numi, freq, temp;
    for (k=0; k<len; k++)
    {
        freq = k * 0.5 / (len - 1);
        zr = cos(-8.0 * atan(1.0) * freq);
        zi = sin(-8.0 * atan(1.0) * freq);
        br = 0.0;
        bi = 0.0;
        for (i=m; i>0; i--)
        {
            re = br;
            im = bi;
            br = (re + b[i]) * zr - im * zi;
            bi = (re + b[i]) * zi + im * zr;
        }
        ar = 0.0;
        ai = 0.0;
        for (i=n; i>0; i--)
        {
            re = ar;
            im = ai;
            ar = (re + a[i]) * zr - im * zi;
            ai = (re + a[i]) * zi + im * zr;
        }
        br = br + b[0];
        ar = ar + 1.0;
        numr = ar * br + ai * bi;
        numi = ar * bi - ai * br;
        den = ar * ar + ai * ai;
        x[k] = numr / den;
        y[k] = numi / den;
        switch (sign)
        {
         case 1:
         {
            temp = sqrt(x[k] * x[k] + y[k] * y[k]);
            y[k] = atan2(y[k], x[k]);
            x[k] = temp;
            break;
         }
         case 2:
         {
            temp = x[k] * x[k] + y[k] * y[k];
            y[k] = atan2(y[k], x[k]);
            x[k] = 10.0 * log10(temp);
         }
        }
    }
}

дҫӢеҰӮпјҡ

ж•°еӯ—зі»з»ҹзҡ„дј йҖ’еҮҪж•°дёәпјҡ
жұӮиҜҘзі»з»ҹзҡ„е№…йў‘е“Қеә”е’Ңзӣёйў‘е“Қеә”пјҢ并画еҮәзӣёеә”зҡ„еӣҫеҪўгҖӮ

жҲ‘жҳҜеҲ©з”ЁQt5.9е®һзҺ°зҡ„гҖӮеҲ—дёҫдё»иҰҒд»Јз ҒгҖӮз»ҳеӣҫеҲ©з”Ёзҡ„жҳҜQCustomPlotгҖӮ

double a[] = {1.0, 0.0, 0.9};
double b[] = {0.0, -0.1};
double x[300];
double y[300];
gain(b, a, 1, 2, x, y, 300, 1);

ж•°еӯ—зі»з»ҹзҡ„е№…йў‘е“Қеә”пјҡ

ж•°еӯ—зі»з»ҹзҡ„зӣёйў‘е“Қеә”пјҡ