changing frequency using fft and ifft not using whole numbers

Question

I know I can change frequency by whole numbers by changing the variable shift but how can I change the frequency using numbers with decimal places like .754 or 1.2345 or 67.456. If I change the variable 'shift' to a non-whole like number like 5.1 I get an error subscript indices must be either positive integers less than 2^31 or logicals from line mag2s = [mag2(shift+1:end), zeros(1,shift)];

Example Code below from question increase / decrease the frequency of a signal using fft and ifft in matlab / octave works with changing the variable shift (but it only works with whole numbers, I need it to work with decimals numbers also).

PS: I'm using octave 3.8.1 which is like matlab and I know I could change the frequency by adjusting the formula in the variable ya but ya will be a signal taken from an audio source (human speech) so it won't be an equation. The equation is just used to keep the example simple. And yes Fs is large due to the fact that signal files used are around 45 seconds long which is why I can't use resample because I get a out of memory error when used.

Here's a animated youtube video example of what I'm trying to get when I use the test equation ya= .5*sin(2*pi*1*t)+.2*cos(2*pi*3*t) and what I'm trying to get happen if I varied the variable shift from (0:0.1:5) youtu.be/pf25Gw6iS1U please keep in mind that ya will be an imported audio signal so I won't have an equation to easily adjust

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

Vav

I'm not sure what you are asking. The frequency shift by shifting this signal is already (probably) a non-integer. The spacing between samples in the fourier domain is 2*Nyquist/N (where N is total number of samples). If you want a closer spacing you can zero pad your input signal. – Ordonez

@Ordonez I know I can change frequency by whole numbers by changing the variable 'shift' but how can I change the frequency using numbers with decimal places like .754 or 1.2345 or 67.456. If I change the variable 'shift' to a non-whole like number 5.1 I get an error subscript indices must be either positive integers less than 2^31 or logicals – Vav

I also added the line where the error comes from – Vav

A general way of frequency shifting is to upsample the signal, mix it with a carrier while suppressing one sideband, then subsample. You can do completely arbitrary frequency shifts that way, the cost is O(N), and you don't deal with FFT :) – Graybeard

@Kuba Ober Thanks, but I'm not exactly sure what you mean or how to program your answer to test this – Vav

@RickT I'm assuming you are targeting a particular frequency? Like you are trying to shift the signal by 5.245 Hz? The FFT is discrete so you can only shift the signal by integer samples, but it is easy to obtain an arbitrary sample spacing by just changing the signal length (while keeping Fs constant). Run your program for t=[0:1/(Fs-1):1] and t=[0:1/(Fs-1):3] and you will see that the amount of frequency change for a five sample shift changes in each case. If this is what you want, I will write a more detailed answer when I get the chance. – Ordonez

@Ordonez I'm trying to frequency shift the imported signal by an exact amount as an example if the imported signal is 3.2hz (please note the imported signals will be human speech audio files) and I want to up shift it up by 5.1hz the new frequency will be 8.3hz. (3.2hz+5.1hz=8.3hz) – Vav

Answer 1

Ok so the question as I understand it is "how do I shift my signal by a specific frequency?"

First let's define Fs which is our sample rate (ie samples per second). We collect a signal which is N samples long. Then the frequency change between samples in the Fourier domain is Fs/N. So taking your example code Fs is 2,000,000 and N is 2,000,000 so the space between each sample is 1Hz and shifting your signal 5 samples shifts it 5Hz.

Now say we want to shift our signal by 5.25Hz instead. Well if our signal was 8,000,000 samples then the spacing would be Fs/N = 0.25Hz and we would shift our signal 11 samples. So how do we get an 8,000,000 sample signal from a 2,000,000 sample signal? Just zero pad it! literally append zeros until it is 8,000,000 samples long. Why does this work? Because you are in essence multiplying your signal by a rectangular window which is equivalent to sinc function convolution in the frequency domain. This is an important point. By appending zeros you are interpolating in the frequency domain (you don't have any more frequency information about the signal you are just interpolating between the previous DTFT points).

We can do this down to any resolution you want, but eventually you'll have to deal with the fact that numbers in digital systems aren't continuous so I recommend just choosing an acceptable tolerance. Lets say we want to be within 0.01 of our desired frequency.

So lets get to actual code. Most of it doesn't change luckily.

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

The final signal is a little over 4Hz which is what we expect. There is some distortion visible from the interpolation, but that should be minimized with a longer signal with a smother frequency domain representation.

Now that I've gone through all of that you may be wondering if there is an easier way. Fortunately for us, there is. We can take advantage of the hilbert transform and fourier transform properties to achieve a frequency shift without ever worrying about Fs or tolerance levels or bin spacing. Namely we know that a time shift leads to a phase shift in the Fourier domain. Well time and frequency are duals so a frequency shift leads to a complex exponential multiplication in the time domain. We don't want to just do a bulk shift of all frequencies because that that will ruin our symmetry in Fourier space leading to a complex time series. So we use the hilbert transform to get the analytic signal which is composed of only the positive frequencies, shift that, and then reconstruct our time series assuming a symmetric Fourier representation.

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

Ordonez

the hilbert transfer works great only if you use single sin waves or single cos waves but if I set FShift =1 and ya = .5*sin(2*pi*1*t)+.2*cos(2*pi*3*t); (I adjusted ya just to try and simulate what the audio import wave may look like) and it seems to fall apart see link to animated youtube video I created showing what happens youtu.be/J5UMrPZvnMA and link to your code I tested it with bit.ly/1lwApNC – Vav

@RickT what are you expecting it to do in that situation? When I looked at it the result with ya= .5*sin(2*pi*1*t)+.2*cos(2*pi*3*t) where FShift = 1, yaShift is exactly yaShift = 0.5*sin(2*pi*2*t)+.2*cos(2*pi*4*t). You have shifted your entire signal by 1hz. – Ordonez

Here's a animated youtube video example of what I'm trying to get when I use the test equation ya= .5*sin(2*pi*1*t)+.2*cos(2*pi*3*t) and what I'm trying to get happen if I varied FShift from (0:0.1:5) youtu.be/pf25Gw6iS1U please keep in mind that ya will be an imported audio signal so I won't have an equation to easily adjust – Vav

Answer 2

You can do this using a fractional delay filter.

First, lets make the code ore workable by letting Matlab handle the conjugate symmetry of the FFT. Just make mag1 and phase1 go to the end . . .

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

Get rid of mag2s and phase2s completely. This simplifies lines 37 and 38 to . .

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

Use the symmetric option of ifft to get Matlb to handle the symmetry for you. You can then drop the forced real, too.

yaifft2 = ifft(yafft2, 'symmetric');         % take inverse fft

With that cleaned up, we can now think of the delay as a filter, e.g.

% ----- changes start here ----- %
shift = 5;
shift_b   = [zeros(1, shift) 1];              % shift amount
shift_a   = 1;

which can be applied as so . . .

mag1s   = filter(shift_b, shift_a, mag1);
phase1s = filter(shift_b, shift_a, phase1);

In this mindset, we can use an allpass filter to make a very simple fractional delay filter

The code above gives the 'M Samples Delay' part of the circuit. You can then add on the fraction using a second cascaded allpass filter . .

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

Winebaum

Thanks so much for your help. You have the variable 'shift' in two places, one as shift = 5 and one as shift = 5.5 is this mislabelled also you have the variables fract_b and fract_a but I can't see where they are used. – Vav

I have been more explicit in the last code block. Look again now – Winebaum

Arbitrarily changing the frequency by using FFT is not generally possible. The size of the FFT will always limit what shifts you can achieve, even if you use fractional delay filters. Try it with a bare sine wave and use an accurate frequency determination algorithm. You'll see that it doesn't work. – Graybeard

@Winebaum The code seems to have issues, I created an animated youtube video of what it does (the blue line is what your code does) at youtu.be/Qir_O8q2bqs and your test code can be found here. bit.ly/1P5zq35 – Vav

Answer 3

Band-limited interpolation using a windowed-Sinc interpolation kernel can be used to change sample rate by arbitrary ratios. Changing the sample rate changes the frequency content of the signal, relative to the sample rate, by the inverse ratio.

Keegan

Thanks, but I'm not exactly sure what you mean or how to program your answer to test this – Vav