MATLAB:如何正确应用 ifft 将 "filtered" 信号带回时域?
MATLAB: How to apply ifft correctly to bring a "filtered" signal back to the time doamin?
我正在尝试获取通过同轴电缆的高斯脉冲输出。我制作了一个代表同轴电缆的矢量;我在网上得到了衰减和相位延迟的信息,并使用欧拉方程创建了一个复杂的阵列。
我对我的高斯向量进行了 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
深入理解FFT的原理对正确使用它很重要。
当你对实信号应用傅里叶变换时,负频率的系数是正频率系数的共轭。当您将 FFT 应用于实数值信号时,您可以从数学上证明应该处于负频率 (-f) 的系数的共轭现在将出现在 (Fsampling-f),其中 Fsampling=1/dt 是采样频率,并且dt 采样周期。这种行为称为混叠,当您将 fft 应用于离散时间信号时会出现,并且采样周期应选择得足够小以使这两个频谱不与香农标准重叠。
当您想对信号应用频率滤波器时,我们说我们保留了频谱的前半部分,因为高频 (>Fsampling/2) 是由于混叠造成的,而不是频谱的特征原始信号。为此,我们在乘以滤波器之前将光谱的后半部分置零。但是,这样做也会丢失原始信号幅度的一半,您将无法使用 ifft 恢复该幅度。选项 'symmetric' 可以通过添加高频 (>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
绘制结果:
如果您还有其他问题,请不要犹豫。祝你好运!
------------编辑---------
离散傅里叶变换的幅度与连续傅里叶变换的幅度不同。如果您有兴趣在频域中显示信号,则需要根据您选择的约定应用归一化。通常,您使用 Dirac delta 函数的傅里叶变换的振幅处处为 1 的约定。
数字 Dirac delta 函数在索引处的振幅为 1,在其他地方为零,导致功率谱在任何地方都等于 1。但是,在您的情况下,时间轴具有采样周期 dt,在这种情况下,数字 Dirac 随时间的积分不是 1,而是 dt。您必须通过将频域信号乘以一个因子 dt(在您的情况下 = 1 皮秒)来规范化频域信号,以遵守约定。您还可以注意到,这使得频域信号与[原始单位乘以时间]同质,这是傅立叶变换的正确单位。
我正在尝试获取通过同轴电缆的高斯脉冲输出。我制作了一个代表同轴电缆的矢量;我在网上得到了衰减和相位延迟的信息,并使用欧拉方程创建了一个复杂的阵列。
我对我的高斯向量进行了 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
深入理解FFT的原理对正确使用它很重要。
当你对实信号应用傅里叶变换时,负频率的系数是正频率系数的共轭。当您将 FFT 应用于实数值信号时,您可以从数学上证明应该处于负频率 (-f) 的系数的共轭现在将出现在 (Fsampling-f),其中 Fsampling=1/dt 是采样频率,并且dt 采样周期。这种行为称为混叠,当您将 fft 应用于离散时间信号时会出现,并且采样周期应选择得足够小以使这两个频谱不与香农标准重叠。
当您想对信号应用频率滤波器时,我们说我们保留了频谱的前半部分,因为高频 (>Fsampling/2) 是由于混叠造成的,而不是频谱的特征原始信号。为此,我们在乘以滤波器之前将光谱的后半部分置零。但是,这样做也会丢失原始信号幅度的一半,您将无法使用 ifft 恢复该幅度。选项 'symmetric' 可以通过添加高频 (>Fsampling/2) 较低系数 ( 我简化了代码以简要解释正在发生的事情,并在第 20 行为您实现了手工制作的对称化。请注意,为了正确显示频谱,我将采样周期从 1 皮秒减少到 100 皮秒: 绘制结果: 如果您还有其他问题,请不要犹豫。祝你好运! ------------编辑--------- 离散傅里叶变换的幅度与连续傅里叶变换的幅度不同。如果您有兴趣在频域中显示信号,则需要根据您选择的约定应用归一化。通常,您使用 Dirac delta 函数的傅里叶变换的振幅处处为 1 的约定。 数字 Dirac delta 函数在索引处的振幅为 1,在其他地方为零,导致功率谱在任何地方都等于 1。但是,在您的情况下,时间轴具有采样周期 dt,在这种情况下,数字 Dirac 随时间的积分不是 1,而是 dt。您必须通过将频域信号乘以一个因子 dt(在您的情况下 = 1 皮秒)来规范化频域信号,以遵守约定。您还可以注意到,这使得频域信号与[原始单位乘以时间]同质,这是傅立叶变换的正确单位。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