изменение частоты с использованием fft и ifft без использования целых чисел

Итак, вопрос, как я понимаю, «как мне сдвинуть мой сигнал на определенную частоту?»

Сначала давайте определим Fs, которая является нашей частотой дискретизации (т.е. выборки в секунду). Мы собираем сигнал длиной N отсчетов. Тогда изменение частоты между отсчетами в области Фурье равно Fs/N. Итак, взяв код вашего примера, Fs равен 2 000 000, а N равен 2 000 000, поэтому расстояние между каждым сэмплом составляет 1 Гц, а сдвиг вашего сигнала на 5 сэмплов сдвигает его на 5 Гц.

Теперь предположим, что вместо этого мы хотим сдвинуть наш сигнал на 5,25 Гц. Что ж, если бы наш сигнал состоял из 8 000 000 сэмплов, тогда интервал был бы Fs/N = 0,25 Гц, и мы бы сдвинули наш сигнал на 11 сэмплов. Итак, как мы можем получить сигнал с 8 000 000 отсчетов из сигнала с 2 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)  ')

Конечный сигнал имеет частоту чуть более 4 Гц, чего мы и ожидали. При интерполяции видны некоторые искажения, но их следует свести к минимуму с помощью более длинного сигнала с более сглаженным представлением в частотной области.

Теперь, когда я прошел через все это, вам может быть интересно, есть ли более простой способ. К счастью для нас, есть. Мы можем воспользоваться преимуществами преобразования Гильберта и преобразования Фурье для достижения сдвига частоты, даже не беспокоясь о 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)  ')

Вы можете сделать это, используя дробный фильтр задержки.

Во-первых, давайте сделаем код работоспособным, позволив Matlab обрабатывать сопряженную симметрию БПФ. Просто сделайте так, чтобы mag1 и phase1 шли в конец. . .

mag1    = mag(2:end);               
phase1  = phase(2:end);     

Полностью избавьтесь от mag2s и phase2s. Это упрощает строки 37 и 38 до . .

magS    = [mag(1)   , mag1s    ];
phaseS  = [phase(1) , phase1s  ];

Используйте опцию symmetric для ifft, чтобы 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);

С таким мышлением мы можем использовать фильтр allpass, чтобы сделать очень простой фильтр с дробной задержкой.

Приведенный выше код дает часть схемы «M Samples Delay». Затем вы можете добавить дробь, используя второй каскадный фильтр allpass. .

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

Интерполяция с ограниченной полосой пропускания с использованием ядра интерполяции с оконным Sinc может использоваться для изменения частоты дискретизации в произвольном соотношении. Изменение частоты дискретизации изменяет частотный состав сигнала относительно частоты дискретизации в обратном отношении.