使用 fft 和 ifft 改变频率而不使用整数

Question

我知道我可以通过更改变量 shift 来更改整数频率，但是如何使用带小数位的数字（如 .754 或1.2345 或 67.456。如果我将变量 'shift' 更改为非整数，如 5.1 我得到一个错误下标索引必须是小于 2^31 的正整数或逻辑行 mag2s = [mag2(shift+1:end), zeros(1,shift)];

下面来自问题 increase / decrease the frequency of a signal using fft and ifft in matlab / octave 的示例代码适用于更改变量 shift（但它只适用于整数，我也需要它适用于小数） .

PS：我使用的是 octave 3.8.1，它类似于 matlab，我知道我可以通过调整变量 ya 中的公式来改变频率，但 ya 会是从音频源（人类语音）获取的信号，因此它不是方程式。该等式仅用于使示例保持简单。是的，Fs 很大，因为使用的信号文件大约 45 秒长，这就是为什么我不能使用 resample，因为我在使用时出现内存不足错误。

这是一个 youtube 动画视频示例，展示了我在使用测试方程 ya= .5*sin(2*pi*1*t)+.2*cos(2*pi*3* 时想要得到的结果t) 以及如果我将变量 shift 从 (0:0.1:5) youtu.be/pf25Gw6iS1U 更改为我想要发生的事情，请记住，你将一个导入的音频信号，所以我没有公式可以轻松调整

clear all,clf

Fs = 2000000;% Sampling frequency
t=linspace(0,1,Fs);

%1a create signal
ya = .5*sin(2*pi*2*t); 

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

mag = abs(ya_fft);
phase = unwrap(angle(ya_fft));
ya_newifft=ifft(mag.*exp(i*phase));

% ----- changes start here ----- %

shift   = 5;                            % shift amount
N       = length(ya_fft);               % number of points in the fft
mag1    = mag(2:N/2+1);                 % get positive freq. magnitude
phase1  = phase(2:N/2+1);               % get positive freq. phases
mag2    = mag(N/2+2:end);               % get negative freq. magnitude
phase2  = phase(N/2+2:end);             % get negative freq. phases

% pad the positive frequency signals with 'shift' zeros on the left
% remove 'shift' components on the right
mag1s   = [zeros(1,shift) , mag1(1:end-shift)];
phase1s = [zeros(1,shift) , phase1(1:end-shift)];

% pad the negative frequency signals with 'shift' zeros on the right
% remove 'shift' components on the left
mag2s   = [mag2(shift+1:end), zeros(1,shift)];
phase2s = [phase2(shift+1:end), zeros(1,shift) ];

% recreate the frequency spectrum after the shift
%           DC      +ve freq.   -ve freq.
magS    = [mag(1)   , mag1s     , mag2s];
phaseS  = [phase(1) , phase1s   , phase2s];


x = magS.*cos(phaseS);                  % change from polar to rectangular
y = magS.*sin(phaseS);
yafft2 = x + i*y;                      % store signal as complex numbers
yaifft2 = real(ifft(yafft2));         % take inverse fft

plot(t,ya,'-r',t,yaifft2,'-b'); % time signal with increased frequency
legend('Original signal (ya)  ','New frequency signal (yaifft2)  ')

Answer 1

您可以使用分数延迟滤波器来做到这一点。

首先，让 Matlab 处理 FFT 的共轭对称性，让代码变得可行。只需让mag1 和phase1 走到最后。 . .

mag1    = mag(2:end);               
phase1  = phase(2:end);     

彻底摆脱mag2s 和phase2s。这将第 37 行和第 38 行简化为 . .

magS    = [mag(1)   , mag1s    ];
phaseS  = [phase(1) , phase1s  ];

使用ifft 的symmetric 选项让Matlb 为您处理对称性。然后你也可以删除强制的real。

yaifft2 = ifft(yafft2, 'symmetric');         % take inverse fft

清理完之后，我们现在可以将延迟视为过滤器，例如

% ----- changes start here ----- %
shift = 5;
shift_b   = [zeros(1, shift) 1];              % shift amount
shift_a   = 1;

可以这样应用。 . .

mag1s   = filter(shift_b, shift_a, mag1);
phase1s = filter(shift_b, shift_a, phase1);

在这种心态下，我们可以使用全通滤波器来制作一个非常简单的分数延迟滤波器

上面的代码给出了电路的“M 个样本延迟”部分。然后，您可以使用第二个级联全通滤波器添加分数。 .

shift = 5.5;
Nw = floor(shift);
shift_b   = [zeros(1, Nw) 1];
shift_a   = 1;

Nf = mod(shift,1);
alpha = -(Nf-1)/(Nf+1);    
fract_b   = [alpha 1];           
fract_a   = [1 alpha];

%// now filter as a cascade . . . 
mag1s   = filter(shift_b, shift_a, mag1);
mag1s   = filter(fract_b, fract_a, mag1s);

Answer 2

好的，据我所知，问题是“如何将信号移动特定频率？”

首先让我们定义 Fs，它是我们的采样率（即每秒采样数）。我们收集一个 N 个样本长的信号。则傅里叶域中样本之间的频率变化为 Fs/N。因此，以您的示例代码为例，Fs 为 2,000,000，N 为 2,000,000，因此每个样本之间的空间为 1Hz，并且将信号移动 5 个样本将其移动 5Hz。

现在假设我们想将信号偏移 5.25Hz。好吧，如果我们的信号是 8,000,000 个样本，那么间距将为 Fs/N = 0.25Hz，我们将把我们的信号移动 11 个样本。那么我们如何从 2,000,000 个样本信号中得到 8,000,000 个样本信号呢？就填零吧！从字面上附加零，直到它有 8,000,000 个样本长。为什么这行得通？因为您本质上是将信号乘以一个矩形窗口，该窗口相当于频域中的 sinc 函数卷积。这是很重要的一点。通过附加零，您将在频域中进行插值（您没有更多关于您只是在先前 DTFT 点之间插值的信号的频率信息）。

我们可以做到这一点，达到您想要的任何分辨率，但最终您将不得不处理数字系统中的数字不连续的事实，因此我建议您选择可接受的容差。假设我们希望在我们想要的频率的 0.01 范围内。

让我们来看看实际的代码。幸运的是，大部分都没有改变。

clear all,clf

Fs = 44100; % lets pick actual audio sampling rate
tolerance = 0.01; % our frequency bin tolerance
minSignalLen = Fs / tolerance; %minimum number of samples for our tolerance

%your code does not like odd length signals so lets make sure we have an 
%even signal length
if(mod(minSignalLen,2) ~=0 )
   minSignalLen = minSignalLen + 1; 
end 

t=linspace(0,1,Fs); %our input signal is 1s long

%1a create 2Hz signal   
ya = .5*sin(2*pi*2*t); 

if (length(ya) < minSignalLen)
    ya = [ya, zeros(1, minSignalLen - length(ya))];
end

df = Fs / length(ya);  %actual frequency domain spacing;

targetFreqShift = 2.32;  %lets shift it 2.32Hz

nSamplesShift = round(targetFreqShift / df);

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

mag = abs(ya_fft);
phase = unwrap(angle(ya_fft));
ya_newifft=ifft(mag.*exp(i*phase));

% ----- changes start here ----- %

shift   = nSamplesShift;                % shift amount
N       = length(ya_fft);               % number of points in the fft
mag1    = mag(2:N/2+1);                 % get positive freq. magnitude
phase1  = phase(2:N/2+1);               % get positive freq. phases
mag2    = mag(N/2+2:end);               % get negative freq. magnitude
phase2  = phase(N/2+2:end);             % get negative freq. phases

% pad the positive frequency signals with 'shift' zeros on the left
% remove 'shift' components on the right
mag1s   = [zeros(1,shift) , mag1(1:end-shift)];
phase1s = [zeros(1,shift) , phase1(1:end-shift)];

% pad the negative frequency signals with 'shift' zeros on the right
% remove 'shift' components on the left
mag2s   = [mag2(shift+1:end), zeros(1,shift)];
phase2s = [phase2(shift+1:end), zeros(1,shift) ];

% recreate the frequency spectrum after the shift
%           DC      +ve freq.   -ve freq. 
magS    = [mag(1)   , mag1s     , mag2s];
phaseS  = [phase(1) , phase1s   , phase2s];


x = magS.*cos(phaseS);                  % change from polar to rectangular
y = magS.*sin(phaseS);
yafft2 = x + i*y;                      % store signal as complex numbers
yaifft2 = real(ifft(yafft2));         % take inverse fft

 %pull out the original 1s of signal
plot(t,ya(1:length(t)),'-r',t,yaifft2(1:length(t)),'-b'); 
legend('Original signal (ya)  ','New frequency signal (yaifft2)  ')

最终信号略高于 4Hz，这是我们所期望的。从插值中可以看到一些失真，但应该使用更长的信号和更平滑的频域表示来最小化失真。

现在我已经完成了所有这些，您可能想知道是否有更简单的方法。对我们来说幸运的是，有。我们可以利用希尔伯特变换和傅里叶变换特性来实现频移，而无需担心 Fs 或容差水平或 bin 间距。也就是说，我们知道时移会导致傅立叶域中的相移。时间和频率是对偶的，因此频移会导致时域中的复指数乘法。我们不想只对所有频率进行批量移动，因为那会破坏我们在傅里叶空间中的对称性，从而导致复杂的时间序列。因此，我们使用希尔伯特变换来获得仅由正频率组成的解析信号，将其移位，然后假设对称傅里叶表示重构我们的时间序列。

Fs = 44100; 
t=linspace(0,1,Fs);
FShift = 2.3 %shift our frequency up by 2.3Hz
%1a create signal
ya = .5*sin(2*pi*2*t);
yaHil = hilbert(ya);  %get the hilbert transform 
yaShiftedHil = yaHil.*exp(1i*2*pi*FShift*t);

yaShifted = real(yaShiftedHil); 
figure
plot(t,ya,'-r',t,yaShifted,'-b') 
legend('Original signal (ya)  ','New frequency signal (yaifft2)  ')

Answer 3

使用窗口 Sinc 插值内核的带限插值可用于以任意比率更改采样率。改变采样率会改变信号的频率成分，相对于采样率，成反比。