C语言实现快速傅里叶变换的算法

「爱情、让人受尽委屈。」 2023-06-08 06:13 8阅读 0赞

快速傅立叶/余弦/正弦变换  
尺寸：一个  
数据长度：2的幂  
抽取：频率  
基数：4，2

### 函数： ###

cdft：复数离散傅立叶变换  
rdft：实离散傅里叶变换  
ddct：离散余弦变换  
ddst：离散正弦变换  
dfct：RDFT（实对称DFT）的余弦变换  
dfst：RDFT（实反对称DFT）的正弦变换

### 函数原型： ###

void cdft(int, int, double *, int *, double *);      // 复数离散傅立叶变换
    void rdft(int, int, double *, int *, double *);      // 实离散傅里叶变换
    void ddct(int, int, double *, int *, double *);      // 离散余弦变换
    void ddst(int, int, double *, int *, double *);      // 离散正弦变换
    void dfct(int, double *, double *, int *, double *); // 实对称DFT
    void dfst(int, double *, double *, int *, double *); // 实反对称DFT

void cdft(int n, int isgn, double *a, int *ip, double *w)
    { 
        
        if (n > (ip[0] << 2)) { 
            makewt(n >> 2, ip, w);
        }
        if (n > 4) { 
            if (isgn >= 0) { 
                bitrv2(n, ip + 2, a);
                cftfsub(n, a, w);
            } else { 
                bitrv2conj(n, ip + 2, a);
                cftbsub(n, a, w);
            }
        } else if (n == 4) { 
            cftfsub(n, a, w);
        }
    }
    
    void rdft(int n, int isgn, double *a, int *ip, double *w)
    { 
        int nw, nc;
        double xi;
        
        nw = ip[0];
        if (n > (nw << 2)) { 
            nw = n >> 2;
            makewt(nw, ip, w);
        }
        nc = ip[1];
        if (n > (nc << 2)) { 
            nc = n >> 2;
            makect(nc, ip, w + nw);
        }
        if (isgn >= 0) { 
            if (n > 4) { 
                bitrv2(n, ip + 2, a);
                cftfsub(n, a, w);
                rftfsub(n, a, nc, w + nw);
            } else if (n == 4) { 
                cftfsub(n, a, w);
            }
            xi = a[0] - a[1];
            a[0] += a[1];
            a[1] = xi;
        } else { 
            a[1] = 0.5 * (a[0] - a[1]);
            a[0] -= a[1];
            if (n > 4) { 
                rftbsub(n, a, nc, w + nw);
                bitrv2(n, ip + 2, a);
                cftbsub(n, a, w);
            } else if (n == 4) { 
                cftfsub(n, a, w);
            }
        }
    }
    
    
    void ddct(int n, int isgn, double *a, int *ip, double *w)
    { 
        int j, nw, nc;
        double xr;
        
        nw = ip[0];
        if (n > (nw << 2)) { 
            nw = n >> 2;
            makewt(nw, ip, w);
        }
        nc = ip[1];
        if (n > nc) { 
            nc = n;
            makect(nc, ip, w + nw);
        }
        if (isgn < 0) { 
            xr = a[n - 1];
            for (j = n - 2; j >= 2; j -= 2) { 
                a[j + 1] = a[j] - a[j - 1];
                a[j] += a[j - 1];
            }
            a[1] = a[0] - xr;
            a[0] += xr;
            if (n > 4) { 
                rftbsub(n, a, nc, w + nw);
                bitrv2(n, ip + 2, a);
                cftbsub(n, a, w);
            } else if (n == 4) { 
                cftfsub(n, a, w);
            }
        }
        dctsub(n, a, nc, w + nw);
        if (isgn >= 0) { 
            if (n > 4) { 
                bitrv2(n, ip + 2, a);
                cftfsub(n, a, w);
                rftfsub(n, a, nc, w + nw);
            } else if (n == 4) { 
                cftfsub(n, a, w);
            }
            xr = a[0] - a[1];
            a[0] += a[1];
            for (j = 2; j < n; j += 2) { 
                a[j - 1] = a[j] - a[j + 1];
                a[j] += a[j + 1];
            }
            a[n - 1] = xr;
        }
    }
    
    
    void ddst(int n, int isgn, double *a, int *ip, double *w)
    { 
        
        int j, nw, nc;
        double xr;
        
        nw = ip[0];
        if (n > (nw << 2)) { 
            nw = n >> 2;
            makewt(nw, ip, w);
        }
        nc = ip[1];
        if (n > nc) { 
            nc = n;
            makect(nc, ip, w + nw);
        }
        if (isgn < 0) { 
            xr = a[n - 1];
            for (j = n - 2; j >= 2; j -= 2) { 
                a[j + 1] = -a[j] - a[j - 1];
                a[j] -= a[j - 1];
            }
            a[1] = a[0] + xr;
            a[0] -= xr;
            if (n > 4) { 
                rftbsub(n, a, nc, w + nw);
                bitrv2(n, ip + 2, a);
                cftbsub(n, a, w);
            } else if (n == 4) { 
                cftfsub(n, a, w);
            }
        }
        dstsub(n, a, nc, w + nw);
        if (isgn >= 0) { 
            if (n > 4) { 
                bitrv2(n, ip + 2, a);
                cftfsub(n, a, w);
                rftfsub(n, a, nc, w + nw);
            } else if (n == 4) { 
                cftfsub(n, a, w);
            }
            xr = a[0] - a[1];
            a[0] += a[1];
            for (j = 2; j < n; j += 2) { 
                a[j - 1] = -a[j] - a[j + 1];
                a[j] -= a[j + 1];
            }
            a[n - 1] = -xr;
        }
    }
    
    
    void dfct(int n, double *a, double *t, int *ip, double *w)
    { 
        int j, k, l, m, mh, nw, nc;
        double xr, xi, yr, yi;
        
        nw = ip[0];
        if (n > (nw << 3)) { 
            nw = n >> 3;
            makewt(nw, ip, w);
        }
        nc = ip[1];
        if (n > (nc << 1)) { 
            nc = n >> 1;
            makect(nc, ip, w + nw);
        }
        m = n >> 1;
        yi = a[m];
        xi = a[0] + a[n];
        a[0] -= a[n];
        t[0] = xi - yi;
        t[m] = xi + yi;
        if (n > 2) { 
            mh = m >> 1;
            for (j = 1; j < mh; j++) { 
                k = m - j;
                xr = a[j] - a[n - j];
                xi = a[j] + a[n - j];
                yr = a[k] - a[n - k];
                yi = a[k] + a[n - k];
                a[j] = xr;
                a[k] = yr;
                t[j] = xi - yi;
                t[k] = xi + yi;
            }
            t[mh] = a[mh] + a[n - mh];
            a[mh] -= a[n - mh];
            dctsub(m, a, nc, w + nw);
            if (m > 4) { 
                bitrv2(m, ip + 2, a);
                cftfsub(m, a, w);
                rftfsub(m, a, nc, w + nw);
            } else if (m == 4) { 
                cftfsub(m, a, w);
            }
            a[n - 1] = a[0] - a[1];
            a[1] = a[0] + a[1];
            for (j = m - 2; j >= 2; j -= 2) { 
                a[2 * j + 1] = a[j] + a[j + 1];
                a[2 * j - 1] = a[j] - a[j + 1];
            }
            l = 2;
            m = mh;
            while (m >= 2) { 
                dctsub(m, t, nc, w + nw);
                if (m > 4) { 
                    bitrv2(m, ip + 2, t);
                    cftfsub(m, t, w);
                    rftfsub(m, t, nc, w + nw);
                } else if (m == 4) { 
                    cftfsub(m, t, w);
                }
                a[n - l] = t[0] - t[1];
                a[l] = t[0] + t[1];
                k = 0;
                for (j = 2; j < m; j += 2) { 
                    k += l << 2;
                    a[k - l] = t[j] - t[j + 1];
                    a[k + l] = t[j] + t[j + 1];
                }
                l <<= 1;
                mh = m >> 1;
                for (j = 0; j < mh; j++) { 
                    k = m - j;
                    t[j] = t[m + k] - t[m + j];
                    t[k] = t[m + k] + t[m + j];
                }
                t[mh] = t[m + mh];
                m = mh;
            }
            a[l] = t[0];
            a[n] = t[2] - t[1];
            a[0] = t[2] + t[1];
        } else { 
            a[1] = a[0];
            a[2] = t[0];
            a[0] = t[1];
        }
    }
    
    
    void dfst(int n, double *a, double *t, int *ip, double *w)
    { 
        int j, k, l, m, mh, nw, nc;
        double xr, xi, yr, yi;
        
        nw = ip[0];
        if (n > (nw << 3)) { 
            nw = n >> 3;
            makewt(nw, ip, w);
        }
        nc = ip[1];
        if (n > (nc << 1)) { 
            nc = n >> 1;
            makect(nc, ip, w + nw);
        }
        if (n > 2) { 
            m = n >> 1;
            mh = m >> 1;
            for (j = 1; j < mh; j++) { 
                k = m - j;
                xr = a[j] + a[n - j];
                xi = a[j] - a[n - j];
                yr = a[k] + a[n - k];
                yi = a[k] - a[n - k];
                a[j] = xr;
                a[k] = yr;
                t[j] = xi + yi;
                t[k] = xi - yi;
            }
            t[0] = a[mh] - a[n - mh];
            a[mh] += a[n - mh];
            a[0] = a[m];
            dstsub(m, a, nc, w + nw);
            if (m > 4) { 
                bitrv2(m, ip + 2, a);
                cftfsub(m, a, w);
                rftfsub(m, a, nc, w + nw);
            } else if (m == 4) { 
                cftfsub(m, a, w);
            }
            a[n - 1] = a[1] - a[0];
            a[1] = a[0] + a[1];
            for (j = m - 2; j >= 2; j -= 2) { 
                a[2 * j + 1] = a[j] - a[j + 1];
                a[2 * j - 1] = -a[j] - a[j + 1];
            }
            l = 2;
            m = mh;
            while (m >= 2) { 
                dstsub(m, t, nc, w + nw);
                if (m > 4) { 
                    bitrv2(m, ip + 2, t);
                    cftfsub(m, t, w);
                    rftfsub(m, t, nc, w + nw);
                } else if (m == 4) { 
                    cftfsub(m, t, w);
                }
                a[n - l] = t[1] - t[0];
                a[l] = t[0] + t[1];
                k = 0;
                for (j = 2; j < m; j += 2) { 
                    k += l << 2;
                    a[k - l] = -t[j] - t[j + 1];
                    a[k + l] = t[j] - t[j + 1];
                }
                l <<= 1;
                mh = m >> 1;
                for (j = 1; j < mh; j++) { 
                    k = m - j;
                    t[j] = t[m + k] + t[m + j];
                    t[k] = t[m + k] - t[m + j];
                }
                t[0] = t[m + mh];
                m = mh;
            }
            a[l] = t[0];
        }
        a[0] = 0;
    }