使用 fft 和 ifft 改变频率不使用整数
changing frequency using fft and ifft not using whole numbers
我知道我可以通过改变变量 shift 来改变频率像 .754 或 1.2345 或 67.456 这样的地方。如果我将变量 'shift' 更改为非整数，如 5.1 我得到一个错误下标索引必须是小于 2^31 的正整数或行 mag2s = [mag2(shift+1:end), zeros(1,shift)];[ 中的 logicals =14=]
下面问题 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 很大，这就是为什么我不能使用重新采样，因为在使用时出现内存不足错误。
这是我在使用测试方程 ya= .5*sin(2*pi*1*t)+.2*cos(2*pi* 3*t) 以及如果我将变量 shift 从 (0:0.1:5) youtu.be/pf25Gw6iS1U 请记住，ya 将是一个导入的音频信号，所以我没有一个可以轻松调整的方程式
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) ')
Band-limited 使用 windowed-Sinc 插值内核的插值可用于按任意比率更改采样率。改变采样率会按反比改变信号的频率内容，相对于采样率。
您可以使用分数延迟滤波器来执行此操作。
首先，让 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 Samples Delay'部分。然后，您可以使用第二个级联全通滤波器添加分数。 .
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);
好的，据我了解，问题是 "how do I shift my signal by a specific frequency?"
首先让我们定义 Fs，这是我们的采样率（即每秒采样数）。我们收集了一个 N 个样本长的信号。那么傅立叶域中样本之间的频率变化为Fs/N。因此，以您的示例代码 Fs 为 2,000,000，N 为 2,000,000，因此每个样本之间的 space 为 1Hz，将信号移动 5 个样本会将其移动 5Hz。
现在假设我们想要将信号偏移 5.25Hz。好吧，如果我们的信号是 8,000,000 个样本，那么间距将为 Fs/N = 0.25Hz，我们会将信号移动 11 个样本。那么我们如何从一个 2,000,000 采样信号中得到一个 8,000,000 采样信号呢？只需零填充即可！从字面上追加零，直到它有 8,000,000 个样本长。为什么这行得通？因为您本质上是将信号乘以矩形 window，这相当于频域中的 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 或容差水平或区间间距。也就是说，我们知道时移会导致傅立叶域中的相移。时间和频率是对偶的，因此频移会导致时域中的复指数乘法。我们不想只对所有频率进行批量移动，因为这会破坏我们在傅立叶 space 中的对称性，从而导致复杂的时间序列。因此，我们使用希尔伯特变换来获取仅由正频率组成的解析信号，对其进行平移，然后假设采用对称傅立叶表示来重建我们的时间序列。
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) ')
我知道我可以通过改变变量 shift 来改变频率像 .754 或 1.2345 或 67.456 这样的地方。如果我将变量 'shift' 更改为非整数，如 5.1 我得到一个错误下标索引必须是小于 2^31 的正整数或行 mag2s = [mag2(shift+1:end), zeros(1,shift)];[ 中的 logicals =14=]
下面问题 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 很大，这就是为什么我不能使用重新采样，因为在使用时出现内存不足错误。
这是我在使用测试方程 ya= .5*sin(2*pi*1*t)+.2*cos(2*pi* 3*t) 以及如果我将变量 shift 从 (0:0.1:5) youtu.be/pf25Gw6iS1U 请记住，ya 将是一个导入的音频信号，所以我没有一个可以轻松调整的方程式
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) ')
Band-limited 使用 windowed-Sinc 插值内核的插值可用于按任意比率更改采样率。改变采样率会按反比改变信号的频率内容，相对于采样率。
您可以使用分数延迟滤波器来执行此操作。
首先，让 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 Samples Delay'部分。然后，您可以使用第二个级联全通滤波器添加分数。 .
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);
好的，据我了解，问题是 "how do I shift my signal by a specific frequency?"
首先让我们定义 Fs，这是我们的采样率（即每秒采样数）。我们收集了一个 N 个样本长的信号。那么傅立叶域中样本之间的频率变化为Fs/N。因此，以您的示例代码 Fs 为 2,000,000，N 为 2,000,000，因此每个样本之间的 space 为 1Hz，将信号移动 5 个样本会将其移动 5Hz。
现在假设我们想要将信号偏移 5.25Hz。好吧，如果我们的信号是 8,000,000 个样本，那么间距将为 Fs/N = 0.25Hz，我们会将信号移动 11 个样本。那么我们如何从一个 2,000,000 采样信号中得到一个 8,000,000 采样信号呢？只需零填充即可！从字面上追加零，直到它有 8,000,000 个样本长。为什么这行得通？因为您本质上是将信号乘以矩形 window，这相当于频域中的 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 或容差水平或区间间距。也就是说，我们知道时移会导致傅立叶域中的相移。时间和频率是对偶的，因此频移会导致时域中的复指数乘法。我们不想只对所有频率进行批量移动，因为这会破坏我们在傅立叶 space 中的对称性，从而导致复杂的时间序列。因此，我们使用希尔伯特变换来获取仅由正频率组成的解析信号，对其进行平移，然后假设采用对称傅立叶表示来重建我们的时间序列。
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) ')