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

Question

我知道我可以通过更改变量 shift 来按整数更改频率，但是如何使用带小数位的数字(如 .754 或)更改频率 1.2345 或 67.456。如果我将变量 'shift' 更改为非整数，如 5.1 我得到一个错误 subscript indices must be either positive integers less than 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 会是从音频源(人类语音)中获取的信号，因此它不会是一个方程式。该等式只是用来保持示例简单。是的，由于使用的信号文件大约有 45 秒长，因此 Fs 很大，这就是为什么我不能使用重新采样，因为在使用时出现内存不足错误。

这是一个动画 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);