基于FFT的频域DSP滤波

n=1024;
n=(1:n)'-1;%'定义时间索引
x=sin(2*pi*1.5*n/n);%输入,每1024点1.5个周期
H=汉宁(129)。*sinc(0.25*(-64:1:64)';%'窗口sinc lpf,fc=pi/4
H=[H./和(H)];%标准化直流增益

y=ifft(fft(x).*fft(h,n));%采样光谱乘积的倒数ft
y=实数(y);%由于数字误差,y有一个很小的虚部
%#根据您的FT/IFT实现,可能需要在这里按n或1/n进行扩展。
情节(Y);
< /代码> 

这是图表:


块开头的故障根本不是我们所期望的。但是,如果考虑到fft(x),它是有意义的。离散傅立叶变换假定信号在变换块内是周期性的。据DFT所知,我们要求对这一时期的一个时期进行转换:


这导致了使用dfts过滤时的第一个重要考虑因素:您实际上正在实现循环卷积,而不是线性卷积。所以当你考虑数学的时候,第一个图表中的“小故障”并不是真正的小故障。那么问题就变成了:有没有办法解决周期性问题?答案是肯定的:使用重叠保存处理。本质上,你计算n长的产品,如上面所述,但只保留n/2分。

nproc=512;
xproc=零(2*nproc,1);%初始化临时缓冲区
idx=1:nproc;%初始化半缓冲区索引
ycorrect=零(2*nproc,1);%初始化目标
对于ctr=1:(长度(x)/nproc)%每次迭代x 512点
xproc(1:nproc)=xproc((nproc+1):结束);%将上一个迭代的下半部分移动到该迭代的上半部分
xproc((nproc+1):end)=x(idx);%用新数据填充此迭代的下半部分
yproc=ifft(fft(xproc)。*fft(h,2*nproc));%计算新缓冲区
ycorrect(idx)=real(yproc((nproc+1):end));%保留新缓冲区的后半部分
idx=idx+nproc;%半步缓冲索引
结束
< /代码> 

以下是ycorrect的图表:


这张图片是有意义的-我们期望从滤波器的启动瞬态,然后结果沉淀到稳态正弦响应。请注意,现在x可以任意长。限制是nproc>2*min(length(x),length(h))。


现在讨论第二个问题:特定的过滤器。在您的循环中,您创建了一个基本上频谱为h=[0(1:511)/512 1(511:-1:1)/512]';if you dohraw=real(ifft(h));plot(hraw),you get:


很难看清,但是在图的最左边缘有一组非零点,然后在极右方的边缘有更多的零点。使用Octave的内置freqzfunction查看我们看到的频率响应(by doingfreqz(hraw)):


幅度响应从高通包络线到零点有很多波动。同样,DFT固有的周期性正在发挥作用。就DFT而言,hrawrepeats一遍又一遍。但是,如果您使用一个周期的hraw,asfreqzdoes,它的光谱与周期性版本有很大的不同。


因此,让我们定义一个新的信号:hrot=[hraw(513:end);hraw(1:512)];we simply rotate the raw dft output to make it continuous within the block.现在,让我们使用freqz(hrot)来查看频率响应。:


好多了。理想的信封就在那里,没有任何涟漪。当然,现在的实现并不那么简单,你必须做一个完全复杂的乘法,而不是仅仅缩放每个复杂的bin,但至少你会得到正确的答案。


请注意,对于速度,您通常会预先计算填充物的DFTh,我将其单独留在循环中,以便更容易与原始数据进行比较。

呃。


滤波传统上是在时域中作为卷积来实现的。你说得对,输入信号和滤波器信号的频谱相乘是等价的。然而,当您使用离散傅立叶变换(DFT)(使用快速傅立叶变换算法实现速度)时,实际上是计算真实频谱的采样版本。这有很多含义,但与滤波最相关的是时域信号是周期性的。

下面是一个例子。考虑一个正弦输入信号x周期内有1.5个周期和一个简单的低通滤波器h.  在matlab/octave语法中:

N = 1024;
n = (1:N)'-1; %'# define the time index
x = sin(2*pi*1.5*n/N); %# input with 1.5 cycles per 1024 points
h = hanning(129) .* sinc(0.25*(-64:1:64)'); %'# windowed sinc LPF, Fc = pi/4
h = [h./sum(h)]; %# normalize DC gain

y = ifft(fft(x) .* fft(h,N)); %# inverse FT of product of sampled spectra
y = real(y); %# due to numerical error, y has a tiny imaginary part
%# Depending on your FT/IFT implementation, might have to scale by N or 1/N here
plot(y);


下面是图表:


块开头的故障根本不是我们所期望的。但是如果你考虑fft(x),这是有道理的。离散傅立叶变换假定信号在变换块内是周期性的。据DFT所知,我们要求对这一时期的一个时期进行转换:


这导致了使用dfts过滤时的第一个重要考虑:您实际上正在实现circular convolution不是线性卷积。所以当你考虑数学的时候,第一个图表中的“小故障”并不是真正的小故障。那么问题就变成了:有没有一种方法可以绕过周期性?答案是肯定的:使用overlap-save processing.  本质上,你计算n长的产品,如上面所述,但只保留n/2分。

Nproc = 512;
xproc = zeros(2*Nproc,1); %# initialize temp buffer
idx = 1:Nproc; %# initialize half-buffer index
ycorrect = zeros(2*Nproc,1); %# initialize destination
for ctr = 1:(length(x)/Nproc) %# iterate over x 512 points at a time
    xproc(1:Nproc) = xproc((Nproc+1):end); %# shift 2nd half of last iteration to 1st half of this iteration
    xproc((Nproc+1):end) = x(idx); %# fill 2nd half of this iteration with new data
    yproc = ifft(fft(xproc) .* fft(h,2*Nproc)); %# calculate new buffer
    ycorrect(idx) = real(yproc((Nproc+1):end)); %# keep 2nd half of new buffer
    idx = idx + Nproc; %# step half-buffer index
end


这是ycorrect:


这张图片是有意义的-我们期望从滤波器的启动瞬态,然后结果沉淀到稳态正弦响应。注意现在X可以任意长。限制是Nproc > 2*min(length(x),length(h)).

现在讨论第二个问题:特定的过滤器。在你的循环中,你创建了一个过滤器,它的光谱本质上是H = [0 (1:511)/512 1 (511:-1:1)/512]';如果你这样做了hraw = real(ifft(H)); plot(hraw)你得到:


很难看清,但是在图的最左边缘有一组非零点,然后在极右方的边缘有更多的零点。使用八度音阶内置freqz函数来查看我们看到的频率响应(通过执行freqz(hraw)):


幅度响应从高通包络线到零点有很多波动。同样,DFT固有的周期性正在发挥作用。就DFT而言,hraw一遍又一遍地重复。但是如果你用一段时间原材料作为数字滤波器频率响应是的,它的光谱和周期版本有很大的不同。

那么让我们定义一个新的信号:hrot = [hraw(513:end) ; hraw(1:512)];我们只需旋转原始DFT输出,使其在块内连续。现在让我们来看看频率响应freqz(hrot):


好多了。理想的信封就在那里,没有任何涟漪。当然,现在实现不是那么简单,你必须做一个完全复杂的乘法fft(hrot)不要只是缩放每个复杂的垃圾箱,但至少你会得到正确的答案。

注意,对于速度,您通常会预先计算填充物的干膜厚度。H我把它单独放在循环中,以便与原来的相比更容易。