MATLAB:如何正确应用 ifft 将 "filtered"信号带回时间域？

Question

我正在尝试获取通过同轴电缆的高斯脉冲输出。我制作了一个代表同轴电缆的矢量；我在网上获得了衰减和相位延迟信息，并使用欧拉方程创建了一个复杂的阵列。

我对我的高斯向量进行了 FFT，并用我的电缆对其进行了卷积。问题是，我不知道如何正确地对卷积进行 iFFT。我在 MathWorks 中阅读了 iFFt，并查看了其他人的问题。有人有类似的问题，在答案中，有人建议删除 n = 2^nextpow2(L) 和 FFT over length(t) 。我能够从中得到更合理的情节，这也解释了为什么会这样。我对是否应该在 iFFt 中使用对称选项感到困惑。这对我的情节产生了很大的影响。我添加对称性的主要原因是因为我在 iFFT 卷积 (timeHF) 中得到了复数。非常感谢您的帮助，谢谢!

clc, clear
Fs = 14E12;          %1 sample per pico seconds
tlim = 4000E-12;
t = -tlim:1/Fs:tlim; %in pico seconds
ag = 0.5;            %peak of guassian 
bg = 0;              %peak location
wg = 50E-12;         %FWHM

x = ag.*exp(-4 .* log(2) .* (t-bg).^2 / (wg).^2); %Gauss. in terms of FWHM
Ly = x;

L = length(t);
%n = 2^nextpow2(L);  %test output in time domain with and without as suggested online
fNum = fft(Ly,L);

frange = Fs/L*(0:(L/2)); %half of the spectrum
fNumMag = abs(fNum/L);   %divide by n to normalize


% COAX modulation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%phase data
mu = 4*pi*1E-7;
sigma_a = 2.9*1E7;
sigma_b = 5.8*1E6;
a = 0.42E-3;
b = 1.75E-3;
er = 1.508;
vf = 0.66;
c = 3E8;
l = 1;
Lso = sqrt(mu) /(4*pi^3/2) * (1/(sqrt(sigma_a)*a) + 1/(b*sqrt(sigma_b)));
Lo = mu/(2*pi) * log(b/a);
%to = l/(vf*c);
to = 12E-9; %measured
phase = -pi*to*(frange + 1/2 * Lso/Lo * sqrt(frange));

%attenuation Data
k1 = 0.34190;
k2 = 0.00377;
len = 1;
mldb = (k1 .* sqrt(frange) + k2 .* frange) ./ 100 .* len ./1E6;
mldb1 = mldb ./ 0.3048; %original eqaution is in inch
tfMag = 10.^(mldb1./-10);
  
% combine to make in complex form
 tfC = [];
     for ii = 1: L/2 + 1
        tfC(ii) = tfMag(ii) * (cosd(phase(ii)) + 1j*sind(phase(ii)));
     end
     
%END ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

%convolute both h and signal
fNum = fNum(1:L/2+1);
convHF = tfC.*fNum;
convHFMag = abs(convHF/L);
timeHF = ifft(convHF, length(t), 'symmetric');  %this is the part im confused about

% Ignore, 
% tfC(numel(fNum)) = 0;
% convHF = tfC.*fNum;
% convHFMag = abs(convHF/n);
% timeHF = ifft(convHF);

%% plotting

% subplot(2, 2, 1);
% plot(t, Ly)
% title('Gaussian input');
% xlabel('time in seconds')
% ylabel('V')
% grid

subplot(2, 2, 1)
plot(frange, abs(tfC(1: L/2 + 1)));
set(gca, 'Xscale', 'log')
title('coax cable model')
xlabel('Hz')
ylabel('|H(s)|V/V')
grid
ylim([0 1.1])

subplot(2, 2, 2);
plot(frange, convHFMag(1:L/2+1), '.-', frange, fNumMag(1:L/2+1)) %make both range and function the same lenght
title('The input signal Vs its convolution with coax');
xlabel('Hz')
ylabel('V')
legend('Convolution','Lorentzian in frequecuency domain');
xlim([0, 5E10])
grid

subplot(2, 2, [3, 4]);
plot(t, Ly, t, timeHF)
% plot(t, real(timeHF(1:length(t)))) %make both range and function the same lenght
legend('Input', 'Output')
title('Signal at the output')
xlabel('time in seconds')
ylabel('V')
grid

Answer 1

深入理解 FFT 的原理以正确使用它很重要。

当您对实信号应用傅里叶变换时，负频率的系数是正频率系数的共轭。当您将 FFT 应用于实数值信号时，您可以从数学上证明应该处于负频率 (-f) 的系数的共轭现在将出现在 (Fsampling-f)，其中 Fsampling=1/dt 是采样频率，并且dt 采样周期。这种行为称为混叠，当您将 fft 应用于离散时间信号时会出现，并且采样周期应选择得足够小以使这两个频谱不与香农标准重叠。

当您想对信号应用频率滤波器时，我们说我们保留了频谱的前半部分，因为高频 (>Fsampling/2) 是由混叠引起的，不是原始信号的特征。为此，我们在乘以滤波器之前将光谱的后半部分置零。但是，这样做也会丢失原始信号幅度的一半，您将无法使用 ifft 恢复该幅度。 “对称”选项可以通过添加高频 (>Fsampling/2) 较低系数 (

我简化了代码以简要解释正在发生的事情，并在第 20 行为您实现了手工制作的对称化。请注意，为了正确显示频谱，我将采样周期从 1 皮秒减少到 100 皮秒:

close all
clc, clear
Fs = 14E10;          %1 sample per pico seconds % CHANGED to 100ps
tlim = 4000E-12;
t = -tlim:1/Fs:tlim; %in pico seconds
ag = 0.5;            %peak of guassian 
bg = 0;              %peak location
wg = 50E-12;         %FWHM
NT = length(t);

x_i = ag.*exp(-4 .* log(2) .* (t-bg).^2 / (wg).^2); %Gauss. in terms of FWHM
fftx_i = fft(x_i);
f = 1/(2*tlim)*(0:NT-1);

fftx_r = fftx_i;
fftx_r(floor(NT/2):end) = 0; % The removal of high frequencies due to aliasing leads to losing half the amplitude
% HER YOU APPLY FILTER
x_r1 = ifft(fftx_r); % without symmetrisation (half the amplitude lost)
x_r2 = ifft(fftx_r, 'symmetric'); % with symmetrisation
x_r3 = ifft(fftx_r+[0, conj(fftx_r(end:-1:2))]); % hand-made symmetrisation

figure();
subplot(211)
hold on
plot(t, x_i, 'r')
plot(t, x_r2, 'r-+')
plot(t, x_r3, 'r-o')
plot(t, x_r1, 'k--')
hold off
legend('Initial', 'Matlab sym', 'Hand made sym', 'No sym')
title('Time signals')
xlabel('time in seconds')
ylabel('V')
grid
subplot(212)
hold on
plot(f, abs(fft(x_i)), 'r')
plot(f, abs(fft(x_r2)), 'r-+')
plot(f, abs(fft(x_r3)), 'r-o')
plot(f, abs(fft(x_r1)), 'k--')
hold off
legend('Initial', 'Matlab sym', 'Hand made sym', 'No sym')
title('Power spectra')
xlabel('frequency in hertz')
ylabel('V')
grid

绘制结果:

如果您还有其他问题，请不要犹豫。祝你好运!

------------ 编辑------------

离散傅立叶变换的振幅与连续傅立叶变换的振幅不同。如果您有兴趣在频域中显示信号，则需要根据您选择的约定应用归一化。通常，您使用狄拉克 delta 函数的傅立叶变换的幅度处处为 1 的约定。

数值 Dirac delta 函数在索引处的振幅为 1，在其他地方为零，导致功率谱在任何地方都等于 1。但是，在您的情况下，时间轴具有采样周期 dt，在这种情况下，数字 Dirac 随时间的积分不是 1，而是 dt。您必须通过将频域信号乘以一个因子 dt(在您的情况下 = 1 皮秒)来规范化频域信号，以遵守约定。您还可以注意到，这使得频域信号与 [原始单位乘以时间] 同质，这是傅里叶变换的正确单位。