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

Question

我知道我可以通过改变变量的整数改变频率偏移,但我怎样才能改变频率使用数字与小数像0.754或1.2345或67.456.如果我将变量'shift'更改为像5.1 这样的非整数,我得到一个错误,下标索引必须是小于2 ^ 31的正整数或来自行mag2s的逻辑= [mag2(shift + 1:end),0 (1,移位)];

示例下面的问题代码在matlab/octave中使用fft和ifft增加/减少信号的频率 与改变变量一起工作(但它只适用于整数,我需要它也可以使用小数数字).

PS:我正在使用octave 3.8.1,就像matlab一样,我知道我可以通过调整变量ya中的公式来改变频率,但是ya将是从音频源(人类语音)中获取的信号,所以它不会是一个等式.该等式仅用于保持示例简单.是的Fs很大,因为使用的信号文件长约45秒,这就是为什么我不能使用重新采样,因为我在使用时出现内存不足错误.

这是一个动画的youtube视频示例,当我使用测试方程ya = .5*sin(2*pi*1*t)+.2*cos(2*pi*3*t)时,我想要得到的内容我想要让发生,如果我改变的变量转变,从(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

好吧，据我了解，问题是“如何将信号转移到特定频率？”

首先，我们定义 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 或容差水平或箱间距。也就是说，我们知道时移会导致傅立叶域中的相移。时间和频率是对偶的，因此频移会导致时域中的复杂指数乘法。我们不想只对所有频率进行批量移位，因为这会破坏傅立叶空间中的对称性，从而导致复杂的时间序列。因此，我们使用希尔伯特变换来获取仅由正频率组成的解析信号，对其进行移位，然后假设对称傅里叶表示来重建我们的时间序列。

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)  ')