二级学院网站建设,简单的景区介绍网页制作代码,沈阳城市建设学院信息与控制工程系,我要设计网计算一个脉冲响应和输入信号的卷积#xff0c;除了使用原始的时域卷积以外#xff0c;还有如下方法#xff1a;
FFT卷积的方法#xff1a;对输入信号#xff08;长度M#xff09;和脉冲响应#xff08;长度N#xff09;分别补零到K#xff08;KMN-1)#xff0c;…计算一个脉冲响应和输入信号的卷积除了使用原始的时域卷积以外还有如下方法
FFT卷积的方法对输入信号长度M和脉冲响应长度N分别补零到KKMN-1)然后分别计算FFT然后相乘最后反FFT取前MN-1个元素作为最终的卷积结果输入信号很长时将输入信号分成一帧一帧使用overlap-add或者overlap-save的方法当脉冲信号和输入信号都很长时可使用均匀分块卷积 均匀分块卷积 均匀分块卷积与频域自适应滤波(FDAF)结合就是WebRTC AEC中线性滤波处理中的算法核心。
在介绍PBFDAF之前我们来看一下均匀分块卷积的计算流程图 我们分几个部分讲解上图的计算流程
1、脉冲响应分块 如上图红色矩形部分将脉冲响应分成P个等长的不重叠的小块每小块的长度为B我们把这些小块叫做子滤波器filter part 1,2...P)将每个小块后面补B个零分别构成2B长度的序列然后进行实数FFT。我们知道实数序列的FFT结果具有对称性因此实数FFT返回B1个点类似numpy的rfft.fft
2、将输入信号分块 如上图红色线框内的部分将输入信号分成不重叠的等长的分块帧分块长度为B通过一个buffer构造重叠长度为B这样输入buffer的长度为2B将输入buffer进行实数FFT得到B1个复数结果。然后将fft结果存入一个长度为P的队列队列进口处存储最新的输入buffer fft结果旧的输入buffer的fft结果从队列的出口扔掉。这个队列有个名字叫Frequency-domain delay line。
3、频域相乘相加和反变换 第三部分如上图红色矩形内将第一部分准备的P个分块脉冲响应的FFT结果分别与第二部分形成的输入buffer fft结果的队列分别相乘然后沿着P的方向求和。得到B1长度的FFT结果然后在进行复数到实数的IFFT如numpy.rfft.ifft)输出是2B个实数序列取后B个元素作为输入block对于的输出。 WebRTC AEC中的分块卷积 % FD block method% ---------------------- Organize dataxk rrin(pos:posN-1);dk ssin(pos:posN-1);xx [xo;xk];xo xk;tmp fft(xx); XX tmp(1:N1);% ---------------------- Filtering XFm(:,1) XX;for mm0:(M-1)mmm1; YFb(:,m) XFm(:,m) .* WFb(:,m);endyfk sum(YFb,2);tmp [yfk ; flipud(conj(yfk(2:N)))];ykt real(ifft(tmp));ykfb ykt(end-N1:end); xk是近端麦克风输入信号要对近端信号进行线性滤波得到估计的回声信号。
xx就是输入bufferxk是输入的N个样本点xo是上一次的输入N个样本点。对输入buffer进行傅里叶变换的到XX将XX存入XFm,XFm就是频域的那个队列
然后将队列中各个输入buffer的fft结果与WFb进行相乘相加。WFb就是存放脉冲响应分块傅里叶变换的结果因为这是自适应滤波所以WFb矩阵的初始值的元素全部是0。M与上文中的P对应N与上文中的B对应。 WebRTC AEC中的PBFDAF
% Partitioned block frequency domain adaptive filtering NLMS and
% standard time-domain sample-based NLMS
%fidfopen(aecFar-samsung.pcm, rb); % Load far end
fidfopen(aecFar.pcm, rb); % Load far end
%fidfopen(farFile, rb); % Load far end
rrinfread(fid,inf,int16);
fclose(fid);
%rrinloadsl(data/far_me2.pcm); % Load far end
%fidfopen(aecNear-samsung.pcm, rb); % Load near end
fidfopen(aecNear.pcm, rb); % Load near end
%fidfopen(nearFile, rb); % Load near end
ssinfread(fid,inf,int16);
%ssin [zeros(1024,1) ; ssin(1:end-1024)];
fclose(fid);
rand(state,13);
fs16000;
multfs/8000;
%rrinrrin(fs*01:round(fs*120));
%ssinssin(fs*01:round(fs*120));% Flags
NLPon0; % NLP
CNon0; % Comfort noise
PLTon1; % Plotting
M 16; % Number of partitions
N 64; % Partition length
L M*N; % Filter length
if fs 8000mufb 0.6;
elsemufb 0.5;
end
%mufb1; alp 0.1; % Power estimation factor alc 0.1; % Coherence estimation factorlenlength(ssin);
wzeros(L,1); % Sample-based TD NLMS
WFbzeros(N1,M); % Block-based FD NLMS
WFbOldzeros(N1,M); % Block-based FD NLMS
YFbzeros(N1,M);
erfbzeros(len,1);zmzeros(N,1);
XFmzeros(N1,M);pn010*ones(N1,1);
pnzeros(N1,1);
NNlen;
Nbfloor(NN/N)-M;start1;
xozeros(N,1);
doxo;
eoxo;for kk1:Nbpos N * (kk-1) start;% FD block method% ---------------------- Organize dataxk rrin(pos:posN-1);dk ssin(pos:posN-1);xx [xo;xk];xo xk;tmp fft(xx); XX tmp(1:N1);% ------------------------ Power estimationpn0 (1 - alp) * pn0 alp * real(XX.* conj(XX));pn pn0;%pn (1 - alp) * pn alp * M * pn0;% ---------------------- Filtering XFm(:,1) XX;for mm0:(M-1)mmm1; YFb(:,m) XFm(:,m) .* WFb(:,m);endyfk sum(YFb,2);tmp [yfk ; flipud(conj(yfk(2:N)))];ykt real(ifft(tmp));ykfb ykt(end-N1:end); % ---------------------- Error estimation ekfb dk - ykfb; erfb(pos:posN-1) ekfb; tmp fft([zm;ekfb]); % FD version for cancelling part (overlap-save)Ek tmp(1:N1);% ------------------------ Adaptation Ek2 Ek ./(M*pn 0.001); % Normalized errorabsEf max(abs(Ek2), threshold);absEf ones(N1,1)*threshold./absEf;Ek2 Ek2.*absEf; % 让EK2限定到thresholdmEk mufb.*Ek2;PP conj(XFm).*(ones(M,1) * mEk); %计算线性相关tmp [PP ; flipud(conj(PP(2:N,:)))];IFPP real(ifft(tmp));PH IFPP(1:N,:);tmp fft([PH;zeros(N,M)]);FPH tmp(1:N1,:);WFb WFb FPH;XFm(:,2:end) XFm(:,1:end-1);endaudiowrite(aec_out.wav, erfb/32768, fs);
