快速卷积与相关
直接计算
利用卷积的定义 \(y[n]=f(n)*g(n)=\sum_{m=0}^{M-1}f[n-m]g[m]\) 若\(f[n],g[n]\)都是实数信号，则需要\(MN\)个乘法
算法：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
/**
* @brief 直接卷积计算
*
* @param x 输入序列
* @param hn 传递函数
* @param m x的长度
* @param n hn的长度
* @param y 输出序列
*/
void direct_conv(double *x, double *hn, int m, int n, double *y)
{
int len1 = m + n - 1;
for (int j = 0; j < len1; j++)
{
y[j] = 0;
for (int i = 0; i < m; i++)
{
if (((j - i) >= 0) && ((j - i) < n))
{
y[j] += x[i] * hn[j - i];
}
}
}
}
快速卷积
用快速傅立叶变换计算线性卷积的算法如下：
1、选择L满足下述条件 \(\begin{cases} L\geq M+N-1 \\ L=2^\gamma &, \text gamma 为正整数 \end{cases}\)
2、将序列 x(n)与 h(n)按如下方式补零,形成长为\(L=2^\gamma\)的序列 \(x(n)=\begin{cases} x(n),& \text n=0,1,...,M-1 \\ 0,& \text n=M,M+1,...,L-1 \end{cases}\)
\[h(n)=\begin{cases} h(n),& \text n=0,1,...,N-1 \\ 0,& \text n=N,N+1,...,L-1 \end{cases}\]3、用FFT算法分别计算x(n)与h(n)的离散傅立叶变换X(k)与H(k) \(X(k)=\sum_{n=0}^{L-1}x(n)e^{-j2\pi nk/L}\)
\[H(k)=\sum_{n=0}^{L-1}h(n)e^{-j2\pi nk/L}\]4、计算 X(k)与H(k)的乘积 \(Y(k)=X(k)H(k)\)
5、用FFT算法计算Y(k)的离散傅立叶反变换，得到卷积y(n） \(y(n)=\frac{1}{L}\sum_{k=0}^{L-1}Y(k)e^{j2\pi nk/L}\)
序列y(n)的前M+N-1点的值就是序列x(n)与y(n)的线性卷积。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
/**
* @brief 使用实数的FFT计算卷积
*
* @param x 开始时存放输入序列x(i)，最后存放线性卷积的值
* @param hn 传递函数hn(i)
* @param m 序列x(i)的长度
* @param n 序列y(i)的长度。
* @param len 线性卷积的长度。len≥m+n-1，且必须是2的整数次幂
*/
void convol(double *x, double *hn, int m, int n, int len)
{
int i, len2;
double t;
for (i = m; i < len; i++)
{
x[i] = 0.0;
}
for (i = n; i < len; i++)
{
hn[i] = 0.0;
}
rfft(x, len);
rfft(hn, len);
len2 = len / 2;
x[0] = x[0] * hn[0];
x[len2] = x[len2] * hn[len2];
for (i = 1; i < len2; i++)
{
t = x[i] * hn[i] - x[len - i] * hn[len - i];
x[len - i] = x[i] * hn[len - i] + x[len - i] * hn[i];
x[i] = t;
}
irfft(x, len);
}
长序列的快速卷积
overlap-save
用重叠保留法(overlap-save)和快速傅立叶变换计算一个长序列和一个短序列的线性卷积。它通常用于数字滤波。
设长序列x(n)的长度为 len，取长度为M的进行分段，序列 h(n)的长度为 N，做L点的FFT，序列 x(n)与 h(n)的卷积定义为：
1、h(n)补偿到L点，后面补零；
2、h(n)做傅里叶变换到H(k)；
3、x(n)分段，重叠前N-1的数据，对于第一段进行前面补零处理；
4、x(n)做傅里叶变换到X(k)；
5、Y(k)=X(k)H(k)；
6、Y(k)做傅里叶反变换到y(k)；
7、保留后L-(N-1)个数据；
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
%overlap save
N = length(h);
L = 64;
M = L-(N-1);
yyy=zeros(1,1024);
nblks=floor(length(x)/M);
hn=[h,zeros(1,L-N)];
H=fft(hn,L);
xx=zeros(1,L);
for j=1:1:nblks
if j==1
xx(1:N-1)=zeros(1,N-1);
xx(N:L)=x(1:M);
else
xx=x((j-1)*M-(N-1)+1:(j-1)*M-(N-1)+L);
end
X=fft(xx,L);
Y=X.*H;
y_=ifft(Y,L);
for i=1:1:M
ii=i+(j-1)*M;
yyy(ii)=y_(i+(N-1));
end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
/**
* @brief 重叠保留法计算卷积
*
* @param x 输入序列
* @param hn 传递函数
* @param y 输出序列
* @param M x的取长度
* @param N hn的长度
* @param L 傅里叶变换长度
* @param len 序列x长度
*/
void overlap_save(double *x, double *hn, double *y, int M, int N, int L, int len)
{
int nblks = len / M;
nblks += len % M ? 1 : 0;
int i, j, i1;
int n2 = L / 2;
double *xx = NULL, *h = NULL, t;
xx = (double *)malloc(L * sizeof(double));
h = (double *)malloc(L * sizeof(double));
memset(y, 0, len * sizeof(double));
memcpy(h, hn, N * sizeof(double));
for (i = N; i < L; i++)
{
h[i] = 0;
}
rfft(h, L);
for (j = 0; j < nblks; j++)
{
if (j == 0)
{
for (i = 0; i < N - 1; i++)
{
xx[i] = 0;
}
for (i = N - 1; i < L; i++)
{
i1 = i - (N - 1);
xx[i] = x[i1];
}
}
else if (j == (nblks - 1))
{
for (i = 0; i < L; i++)
{
i1 = j * M - (N - 1) + i;
xx[i] = i1 < len ? x[i1] : 0;
// xx[i]=0;
}
}
else
{
for (i = 0; i < L; i++)
{
i1 = j * M - (N - 1) + i;
xx[i] = x[i1];
}
}
rfft(xx, L);
/*复数乘法 r*h */
xx[0] = xx[0] * h[0];
xx[n2] = xx[n2] * h[n2];
for (i = 1; i < n2; i++)
{
t = h[i] * xx[i] - h[L - i] * xx[L - i];
xx[L - i] = h[i] * xx[L - i] + h[L - i] * xx[i];
xx[i] = t;
}
irfft(xx, L);
/*保存结果数据*/
if (j == (nblks - 1))
{
for (i = 0; i < M; i++)
{
i1 = i + j * M;
if (i1 < len)
{
y[i1] = xx[i + (N - 1)];
}
}
}
else
{
for (i = 0; i < M; i++)
{
i1 = i + j * M;
y[i1] = xx[i + (N - 1)];
}
}
}
free(xx);
free(h);
}
替换为risc-v的fft库，去除内存分配
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
/**
* @brief 重叠保留法计算卷积
*
* @param x 输入序列
* @param hn 传递函数
* @param y 输出序列
* @param M x的取长度
* @param N hn的长度
* @param L 傅里叶变换长度
* @param len 序列x长度
*/
void overlap_save2(float *x, float *hn, float *y, int M, int N, int L, int len)
{
int nblks = len / M;
nblks += len % M ? 1 : 0;
int i, j, i1;
int n2 = L / 2;
float *xx = NULL, *h = NULL, t;
float *xx_out=NULL, *h_out=NULL;
// xx = (float *)malloc(L * sizeof(float));
// h = (float *)malloc(L * sizeof(float));
xx = temp1;
h = temp2;
xx_out = temp3;
h_out = temp4;
memset(y, 0, len * sizeof(float));
memcpy(h, hn, N * sizeof(float));
for (i = N; i < L; i++)
{
h[i] = 0;
}
// rfft(h, L);
riscv_rfft_fast_f32(&S, h, h_out, 0);
for (j = 0; j < nblks; j++)
{
if (j == 0)
{
for (i = 0; i < N - 1; i++)
{
xx[i] = 0;
}
for (i = N - 1; i < L; i++)
{
i1 = i - (N - 1);
xx[i] = x[i1];
}
}
else if (j == (nblks - 1))
{
for (i = 0; i < L; i++)
{
i1 = j * M - (N - 1) + i;
xx[i] = i1 < len ? x[i1] : 0;
}
}
else
{
for (i = 0; i < L; i++)
{
i1 = j * M - (N - 1) + i;
xx[i] = x[i1];
}
}
// rfft(xx, L);
riscv_rfft_fast_f32(&S, xx, xx_out, 0);
/*复数乘法 x*h */
xx_out[0] = xx_out[0] * h_out[0];
xx_out[1] = xx_out[1] * h_out[1];
for (i = 1; i < n2; i++)
{
t = h_out[2 * i] * xx_out[2 * i] - h_out[2 * i + 1] * xx_out[2 * i + 1];
xx_out[2 * i + 1] = h_out[2 * i] * xx_out[2 * i + 1] + h_out[2 * i + 1] * xx_out[2 * i];
xx_out[2 * i] = t;
}
// irfft(xx, L);
riscv_rfft_fast_f32(&S, xx_out, xx, 1);
/*保存结果数据*/
if (j == (nblks - 1))
{
for (i = 0; i < M; i++)
{
i1 = i + j * M;
if (i1 < len)
{
y[i1] = xx[i + (N - 1)];
}
}
}
else
{
for (i = 0; i < M; i++)
{
i1 = i + j * M;
y[i1] = xx[i + (N - 1)];
}
}
}
// free(xx);
// free(h);
}
这里只有切换到快速的FFT计算之后才会比直接计算快很多。
overlap-add
设长序列x(n)的长度为 len，取长度为M的进行分段，序列 h(n)的长度为 N，做L点的FFT，序列 x(n)与 h(n)的卷积定义为：
1、h(n)补偿到L点，后面补零；
2、h(n)做傅里叶变换到H(k)；
3、x(n)分段，取前M个数据，后L-M个补零
4、x(n)做傅里叶变换到X(k)；
5、Y(k)=X(k)H(k)；
6、Y(k)做傅里叶反变换到y(k)；
7、对于混叠的(N-1)点做相加处理；
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
%overlap add
N = length(h);
M = 40;
L = 64;
yy=zeros(1,1024);
nblks=floor(length(x)/M);
hn=[h,zeros(1,L-N)];
H=fft(hn,L);
for j=1:1:nblks
xx=[x((j-1)*M+1:j*M),zeros(1,L-M)];
X=fft(xx,L);
Y=X.*H;
y_=ifft(Y,L);
for i=1:1:N+M-1
ii=i+(j-1)*M;
yy(ii)=yy(ii)+y_(i);
end
end
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
/**
* @brief 重叠相加法计算卷积
*
* @param x 输入序列
* @param hn 传递函数
* @param y 输出序列
* @param M x的取长度
* @param N hn的长度
* @param L 傅里叶变换长度
* @param len 序列x长度
*/
void overlap_add(double *x, double *hn, double *y, int M, int N, int L, int len)
{
int nblks = len / M;
nblks += len % M ? 1 : 0;
int i, j, i1;
int n2 = L / 2;
double *xx = NULL, *h = NULL, t;
xx = (double *)malloc(L * sizeof(double));
h = (double *)malloc(L * sizeof(double));
memset(y, 0, len * sizeof(double));
memcpy(h, hn, N * sizeof(double));
for (i = N; i < L; i++)
{
h[i] = 0;
}
rfft(h, L);
for (j = 0; j < nblks; j++)
{
if (j == (nblks - 1))
{
for (i = 0; i < M; i++)
{
i1 = i + j * M;
xx[i] = i1 < len ? x[i1] : 0;
}
}
else
{
for (i = 0; i < M; i++)
{
i1 = i + j * M;
xx[i] = x[i1];
}
}
for (i = M; i < L; i++)
{
xx[i] = 0;
}
rfft(xx, L);
/*复数乘法 */
xx[0] = xx[0] * h[0];
xx[n2] = xx[n2] * h[n2];
for (i = 1; i < n2; i++)
{
t = h[i] * xx[i] - h[L - i] * xx[L - i];
xx[L - i] = h[i] * xx[L - i] + h[L - i] * xx[i];
xx[i] = t;
}
irfft(xx, L);
/*保存结果数据*/
if (j == (nblks - 1))
{
for (i = 0; i < M + (N - 1); i++)
{
i1 = i + j * M;
if (i1 < len)
{
y[i1] += xx[i];
}
}
}
else
{
for (i = 0; i < M + (N - 1); i++)
{
i1 = i + j * M;
y[i1] += xx[i];
}
}
}
free(xx);
free(h);
}
注释：
这里的复数乘法为rfft输出的格式影响，如果复数为：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
// 复数乘法
Complex complexMultiply(Complex *a, Complex *b)
{
Complex result;
result.real = a->real * b->real - a->imag * b->imag;
result.imag = a->real * b->imag + a->imag * b->real;
return result;
}
// 复数除法
Complex complexDivide(Complex *a, Complex *b)
{
Complex result;
double denominator = b->real * b->real + b->imag * b->imag;
if (denominator == 0.0)
{
// fprintf(stderr, "Error: Division by zero.\n");
result.real = 0.0;
result.imag = 0.0;
}
else
{
result.real = (a->real * b->real + a->imag * b->imag) / denominator;
result.imag = (a->imag * b->real - a->real * b->imag) / denominator;
}
return result;
}
应用时候
以上三种方法皆可用来计算卷积，其差别在于所需总体乘法量不同。基于运算量以及效率的考量，在计算卷积时，通常会选择所需总体乘法量较少的方法。
以下根据f(n)和h(n)长度(M,N)分成5类，输入序列长M，传递函数长N，并列出适合使用的方法：
1、M，N很小，直接计算
2、M，N差不多，快速傅里叶计算
3、M，N相差很大，采用分段卷积
在实际应用中，这里不能使用rfft，计算太慢了，使用arm库开源的查表fft还行。
快速相关
使用快速傅里叶变换计算两个有限长序列的线性相关。
1、选择长度L
L >= M+N-1 , L=2^n
2、将序列x(n)和y(n)按如下方式补零。
x长度为M，往后补零到L
y长度为N，前M-1补零，中间填充，M+N-1~L也补零
3、计算L点长的FFT变换
x –> X y –> Y
4、计算X’ (共轭复数) 与Y的乘积
5、傅里叶反变换
6、前M+N-1点就是x与y的线性相关。
本文由作者按照 CC BY 4.0 进行授权