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

Question

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

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

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

请参阅下面的简单示例:我将应用 FFT 的信号是大约 60 秒长、44100 个样本(大型数组)的人类音频信号，因此我正在尝试查看是否可以使用/更改 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.

Answer 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);

我编写了这个函数，使其使用代表自波的单个周期的采样信号。

以下函数使用 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中使用了余弦函数。它当然应该是一个复指数。但我只是采取捷径假设被逆变换的数据是共轭对称的。如果正变换的输入是实数，情况总是如此(逆变换的输出也必须是实数，由于共轭对称性，两个频率的复分量相互抵消)。