更改快速傅里叶逆变换 (ifft) 以使用任意波形而不是正弦波来创建新信号

Changing the inverse fast Fourier transform (ifft) to use an arbitrary waveform instead of sine waves to create a new signal

我知道快速傅里叶逆变换 (ifft) 将从通过对信号执行 fft 获得的数据中的多个正弦波求和。 有没有办法使用新型的快速傅立叶逆变换 (ifft) 使用任意波形而不是仅使用正弦波来创建信号?

我不是要重新创建原始信号。我正在尝试使用一种新型的快速傅里叶逆变换 (ifft) 创建一个新信号,该变换使用给定的任意波形,该波形基于从源信号的 fft 计算得出的(频率、幅度、相位)数据。

任意波形是一种采样信号,它将取代 fft 中使用的正弦波的一个周期。也就是说,信号要根据fft给出的值进行缩放、重复和移位。

请参阅下面的简单示例:我将应用 FFT 的信号是 44100 个样本(大型阵列)时长约 60 秒的人类音频信号,因此 我想看看我是否可以使用 /以某种方式改变 ifft 命令以使用/基于任意波形创建新信号。

PS: 我正在使用类似于 Matlab 的 Octave 4.0,用于创建新信号的任意波形信号将被更改以创建不同的信号。

clear all,clf reset, clc,tic

fs=44100 % Sampling frequency
len_of_sig=2; %length of signal in seconds
t=linspace(0,2*pi*len_of_sig,fs*len_of_sig);
afp=[.5,2.43,pi/9;.3,3,pi/2;.3,4.3,pi/3];  %represents Amplitude,frequency,phase data array

%1 create source signal
ya=0;
for zz=1:size(afp,1)
  ya = ya+afp(zz,1)*sin(afp(zz,2)*t+afp(zz,3));
end

%2 create source frequency domain data
ya_fft = fft(ya);

%3 rebuild original source signal
mag = abs(ya_fft);
phase = unwrap(angle(ya_fft));
ya_newifft=ifft(mag.*exp(i*phase));
ifft_sig_combined_L1=ifft(mag.*exp(i*phase),length(ya_newifft)); 

%4 %%%-----begin create arbitrary waveform to use ---- 
gauss = @(t, t0, g) exp(-((t-t0)/g).^2); % a simple gaussian

t_arbitrary=0:1:44100; % sampling
t_arbitrary_1 = 10000; % pulses peak positions (s)
t_arbitrary_2 = 30000; % pulses peak positions (s)

g = 2000; % pulses width (at 1/e^2) (s)

lilly = gauss(t_arbitrary, t_arbitrary_1, g) - (.57*gauss(t_arbitrary, t_arbitrary_2, g)); %different amplitude peaks
%%%%-----End arbitrary waveform to use---- 

%5 plot
t_sec=t./(2*pi); %converts time in radians to seconds
t_arbitrary_sec=t_arbitrary./length(lilly); %converts time in radians to seconds

subplot(4,1,1);
plot(t_sec,ya,'r')
title('1) source signal')

subplot(4,1,2);
plot(t_sec,ifft_sig_combined_L1)
title('2) rebuilt source signal using ifft')

subplot(4,1,3);
plot(t_arbitrary_sec,lilly,'r')
title('3) arbitrary waveform used to create new signal')

在下面添加了一个带有简单信号的工作流程图,看看是否能更好地解释它:

Section 1) The audio signal is read into an array
Section 2) FFT is done on the signal
Section 3 Red) Normally Inverse FFT uses sin waves to rebuild the signal see signal in red
Section 3 Blue) I want to use an arbitrary signal wave instead to rebuild the signal using the FFT data calculated in (Section 2)
Section 4) New signals created using a new type of Inverse FFT (Section 3).
Please note the new type of Inverse FFT final signal (in blue ) must use the FFT data taken from the original signal.
The signal Sample rate tested should be 44100 and the length of the signal in seconds should be 57.3 seconds long.  I use these numbers to test that the array can handle large amounts and that the code can handle non even numbers in seconds.

IFFT 只是实现 IDFT 的一种方式。 IDFT 只是孔径为整数周期的正弦波形的加权和。

如果需要,您可以采用几乎任何 DFT 或 IDFT 算法或源代码,并将 sin() 函数替换为您想要用于波形合成的任何其他函数。如果您愿意,您甚至可以为不同的频率使用不同的波形,或者将合成频率更改为非整数周期孔径。

(逆)快速傅立叶变换依赖于正弦函数的特殊属性,它允许以低得多的计算成本在时域和频域之间移动 (O(n.log(n)) ) 而不是任意波形 (O(n^2))。如果将基本波形从正弦曲线更改为其他波形,通常,您将无法再获得 FFT 的计算优势。

在您的情况下,听起来您可能想要生成一个与原始信号具有相同频谱但不一定在时域具有相同包络的信号。正如我认为您已经在代码中实现的那样,最简单的方法就是简单地更改 FFT 输出中每个频率仓的相位,然后进行逆 FFT。这将生成一个时域信号,其外观与您的输入信号截然不同,但每个频率仓中的功率相同。

您可能需要牢记的一个微妙之处是如何更改相位以使输出信号保持实值,而不是涉及复数。您可以将傅里叶变换视为产生一组正频率和负频率。对于实值信号,相应的正频率和负频率必须具有彼此复共轭的幅度。假设您的输入信号是实值,您的 FFT 将已经具有此 属性,但您需要安排您应用的随机相位仍然保持正负频率之间的这种关系。实际上,这意味着您只有大约一半的随机相位可供选择 - none 用于零频率仓,每个正频率(第一个 n/2 条目) FFT,其他阶段使得入口 k 处的相位是 (n-k) 处入口的 -1 倍。

让我们从函数 lilly 开始,该函数采用频率、振幅和相位(所有标量)以及信号长度 N,并按预期计算正弦波逆 DFT(见下面的注释 2):

function out = lilly(N,periods,amp,phase)
persistent t
persistent oneperiod
if numel(t)~=N
   disp('recomputung "oneperiod"');
   t = 0:N-1;
   oneperiod = cos(t * 2 * pi / N);
end
p = round(t * periods + phase/(2*pi)*N);
p = mod(p,N) + 1;
out = amp * oneperiod(p);

我编写了这个函数,它使用代表 since 波的单个周期的采样信号。

以下函数使用 lilly 函数计算逆 DFT(见下面的注释 1):

function out = pseudoifft(ft)
N = length(ft);
half = ceil((N+1)/2);
out = abs(ft(1)) + abs(ft(half)) * ones(1,N);
for k=2:half-1
   out = out + lilly(N,k-1,2*abs(ft(k)),angle(ft(k)));
end
out = out/N;

现在我测试以验证它确实计算了逆 DFT:

>> a=randn(1,256);
>> b=fft(a);
>> c=pseudoifft(b);
recomputung "oneperiod"
>> max(abs(a-c))
ans =  0.059656
>> subplot(2,1,1);plot(a)
>> subplot(2,1,2);plot(c)

误差比较大,由于round函数:我们对信号进行子采样而不是插值。如果您需要更高的精度(我认为不太可能),您应该使用 interp1 而不是使用 round(p).

进行索引

接下来,我们用您的示例信号替换 lilly 函数中的正弦:

function out = lilly(N,periods,amp,phase)
persistent t
persistent oneperiod
if numel(t)~=N
   disp('recomputung "oneperiod"');
   t = 0:N-1;
   %oneperiod = cos(t * 2 * pi / N);
   gauss = @(t,t0,g) exp(-((t-t0)/g).^2); % a simple gaussian
   t1 = N/4;   % pulses peak positions (s)
   t2 = 3*N/4; % pulses peak positions (s)
   g = N/20;   % pulses width (at 1/e^2) (s)
   oneperiod = gauss(t,t1,g) - (.57*gauss(t,t2,g)); %different amplitude peaks
   oneperiod = circshift(oneperiod,[1,-round(N/4)]); % this will make it look more like cos
end
p = round(t * periods + phase/(2*pi)*N);
p = mod(p,N) + 1;
out = amp * oneperiod(p);

函数 pseudoifft 现在创建一个由您的基础组成的函数:

>> c=pseudoifft(b);
recomputung "oneperiod"
>> subplot(2,1,2);plot(c)

让我们看一个更简单的输入:

>> z=zeros(size(a));
>> z(10)=1;
>> subplot(2,1,1);plot(pseudoifft(z))
>> z(19)=0.2;
>> subplot(2,1,2);plot(pseudoifft(z))


注1:在你的问题中你特别要求使用FFT。 FFT 只是一种计算正向和反向 DFT 的有效方法。上面的代码在 O(n^2) 中计算逆 DFT,FFT 将在 O(n log n) 中计算相同的结果。不幸的是,FFT 是一种基于 DFT 中使用的复指数属性构建的算法,如果用任何其他函数替换该复指数,则无法使用相同的算法。

注2:我在反DFT中使用余弦函数。它当然应该是一个复杂的指数。但我只是走捷径,假设被逆变换的数据是共轭对称的。如果正向变换的输入是实数(逆变换的输出也必须是实数,两个频率的复分量由于共轭对称性而抵消),情况总是如此。