时间卷积与频率乘法

Question

我正在尝试获取某个频段的功率，但我想在时域而不是频域中执行此操作。问题 - 波段非常紧，因此使用简单的过滤器会产生重叠的“尾巴”。

让 [a1 a2]Hz 是我要计算功率的频段。我可以想象我将频域乘以一个矩形信号，因此我可以及时得到它，所以我可以及时进行卷积。

代码是：（matlab） x - 时间信号，X=fft(x)，W - 频率窗口，w=ifft(W)

filteredX=X.*W;
Fx=ifft(filteredX);
Fx2=conv(x,w,'same');

信号的结果不同。 虽然经过过滤的X 显示了正确的频谱，但 Fx2 的 fft（卷积的结果）却完全不同。

有什么建议吗？

编辑：

根据 EitanT 的建议（谢谢），我使用了以下代码：

im = fix(255 * rand(500,1)); 
mask = ones(4,1) / 16; 

% # Circular convolution
resConv = conv(im, mask);

% # Discrete Fourier transform
M = size(im, 1) + size(mask, 1);
resIFFT = ifft(fft(im, M) .* fft(mask, M));
% # Not needed any more - resIFFT = resIFFT(1:end-1);  % # Adjust dimensions

% # Check the difference
max(abs(resConv(:) - resIFFT(:)))

效果很好，但我不能使用它，所以我不得不更改与尺寸问题有关的部分并得到以下（请参阅 cmets）：

 im = fix(255 * rand(500,1)); 
mask = ones(4,1) / 16; 

% # Circular convolution
resConv = conv(im, mask,'same'); % # instead of conv(im, mask)

% # Discrete Fourier transform
M = size(im, 1) % # Instead of: M = size(im, 1) + size(mask, 1);
resIFFT = ifft(fft(im, M) .* fft(mask, M));
resIFFT = resIFFT(1:end-1);  % # Adjust dimensions

% # Check the difference
max(abs(resConv(:) - resIFFT(:)))

虽然我希望得到相同的结果，但现在差异要大得多。

Answer 1

这是一个很好的例子： https://www.youtube.com/watch?v=iUafo2UZowE

编辑：

只是不要忘记填充。对于一维信号，它是这样的：

lh=length(h);
lx=length(x);

h=[h zeros(1,lx-1)];
x=[x zeros(1,lh-1)];

剩下的很简单：

H=fft(h);
X=fft(x);

Y=H.*X;

y=ifft(Y);

stem(y)

如果您想绘制响应与时间的关系：

假设nh和nx是h和x样本的对应时间，那么响应样本的对应时间可以这样计算：

n=min(nh)+min(nx):+max(nh)+max(nx);

最后

stem(n,y)

Answer 2

Eitan 早期的verification of the convolution 定理非常出色。在这里，我想演示用于过滤特定频段的时域卷积，并展示它等效于频域乘法。

除了知道如何运行fft 和ifft 函数之外，还需要两个部分来完成这项工作：

1. 确定fft 函数返回的每个傅立叶分量的频率 (Hz)
2. 在fft 之前填充信号，在ifft 或conv 函数之后对其进行裁剪。

下面的代码生成一个随机信号，同时执行时域和频域滤波，并在图中显示了等价性：

Nsamp = 50;  %number of samples from signal X
X = randn(Nsamp,1);  %random Gaussian signal
fs = 1000;     %sampling frequency

NFFT = 2*Nsamp;  %number of samples to be used in fft (this will lead to padding)

bandlow = 250;  %lower end of filter band (just an example range)
bandhigh = 450; %upper end of filter band

%construct a frequency axis so we know where Fourier components are

if mod(NFFT,2) == 0 %if there is an even number of points in the FFT
    iNyq = NFFT/2;     %index to the highest frequency
    posfqs = fs*(0:iNyq)/NFFT;  %positive frequencies
    negfqs = -posfqs(end-1:-1:2);  %negative frequencies

else  %if there is an odd number of point in the FFT
    iNyq = floor(NFFT/2);     %index to the highest frequency
    posfqs = fs*(0:iNyq)/NFFT;  %positive frequencies
    negfqs = -posfqs(end:-1:2);  %negative frequencies
end

fqs = [posfqs'; negfqs'];   %concatenate the positive and negative freqs

fftX = fft(X,NFFT);  % compute the NFFT-point discrete Fourier transform of X
                     % becuse NFFT > Nsamp, X is zero-padded

% construct frequency-space mask for the desired frequency range
W = ones(size(fftX)) .* ( and ( abs(fqs) >=  bandlow, abs(fqs) < bandhigh) );

fftX_filtered = fftX.*W; %multiplication in frequency space
X_mult_filtered = ifft( fftX_filtered, NFFT);  %convert the multiplicatively-filtered signal to time domain

w = ifft(W,NFFT); %convert the filter to the time domain
X_conv_filtered = conv(X, w); %convolve with filter in time domain

%now plot and compare the frequency and time domain results
figure; set(gcf, 'Units', 'normalized', 'Position', [0.45 0.15 0.4 0.7])

subplot(3,1,1); 
plot(X)
xlabel('Time Samples'); ylabel('Values of X'); title('Original Signal')

subplot(3,1,2); 
plot(X_conv_filtered(1:Nsamp))
xlabel('Time Samples'); ylabel('Values of Filtered X'); title('Filtered by Convolution in Time Domain')

subplot(3,1,3); 
plot(X_mult_filtered(1:Nsamp))
xlabel('Time Samples'); ylabel('Values of Filtered X'); title('Filtered by Multiplication in Frequency Domain')

for sp = 1:3; subplot(3,1,sp); axis([0 Nsamp -3 3]); end