使用FFT和IFFT改变频率不使用整数

Question

我知道我可以通过改变可变转变的整数改变频率，但我怎样才能改变频率使用数字与小数像0.754或1.2345或67.456。如果我改变可变“转移”到非整像数等5.1我得到一个错误标索引必须是正整数小于2 ^从管线mag2s = [MAG2 31或逻辑值（移+1：结束），零（1，移位）];从问题increase/decrease the frequency of a signal using fft and ifft in matlab/octave不断变化的变量转变作品（只是到整数的作品，我需要它带有小数数字也工作）以下

示例代码。

PS：我使用八度3.8.1这就好比是MATLAB和我知道我可以通过在变量调整式改变频率雅但亚将来自音频源取的信号（人类的言语），所以它不会是一个等式。该等式仅用于保持示例简单。而且Fs很大，因为使用的信号文件大约45秒长，这就是为什么我不能使用resample，因为使用时出现内存不足错误。

这里是一个动画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个样本。那么我们如何从200万个采样信号中获得800万个采样信号呢？只需零填充它！从字面上追加零直到800万个样本。为什么这个工作？因为你本质上是将信号乘以一个等效于频域中的辛函数卷积的矩形窗口。这是重要的一点。通过附加您在频域上插零（你没有对信号你只是以前的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) ') 

的最终信号略微超过4赫兹这是我们所期望的。从插值中可以看到一些失真，但是应该用具有扼杀频域表示的较长信号来使其最小化。

现在，我已经通过了这一切以后，你可能会想，如果有一个更简单的方法。对我们来说幸运的是，有。我们可以利用希尔伯特变换和傅里叶变换特性来实现频移，而不必担心F或公差水平或箱间距。也就是说，我们知道时间偏移会导致傅立叶域的相移。好的时间和频率是双重的，所以频移会导致时域中复杂的指数倍增。我们不想仅仅做所有频率的整体移位，因为这会破坏傅立叶空间中的对称性，导致复杂的时间序列。因此，我们使用希尔伯特变换来获得仅由正频率组成的分析信号，然后对其进行偏移，然后重构假设对称傅立叶表示的时间序列。

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 2

带限内插使用加窗非晶硅nc插值内核可用于以任意比率改变采样率。改变采样率可以改变信号的频率成分相对于采样率的倒数比例。

Answer 3

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

首先，通过让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);