使用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);