Different Fourier transform (fft) of the same signal with different sampling frequency

The fft function in Matlab returns different output for different sampling frequency of the same signal.

Shouldn’t the Fourier transform be insensitive to different sampling frequency of the same signal?

clear all;
M = 0.1;
c=0.15;

% Signal 1
fs1 = 50; %sampling frequency
rs1 = 1/fs1;  
r1 = 0:rs1:2-rs1; 
L1 = length(r1); %signal length
% Gaussian disribution 
M1 = M * exp(-1/2*(r1/c).^2);
n1 = 2^nextpow2(L1);
%Convert the Gaussian pulse to the frequency domain.
FFT1 = fft(M1); 
%Define the frequency domain 
f1 = 0:(fs1/n1):(fs1/2-fs1/n1);

% Signal 2
fs2 = 500; %sampling frequency
rs2 = 1/fs2; 
r2 = 0:rs2:2-rs2; 
L2 = length(r2); %signal length
% Gaussian disribution
M2 = M * exp(-1/2*(r2/c).^2);
n2 = 2^nextpow2(L2);
%Convert the Gaussian pulse to the frequency domain.
FFT2 = fft(M2); 
%Define the frequency domain 
f2 = 0:(fs2/n2):(fs2/2-fs2/n2);

%%
figure;
plot(r2,M2,'.','color','red'); hold on;
plot(r1,M1,'x','color','blue'); hold on;
xlabel('r');
ylabel('M');
title('Signal');
legend('Sampling frequency 500','Sampling frequency 50');

%%
figure;
plot(f2,abs(FFT2(1:n2/2)),'color','red'); hold on;
plot(f1,abs(FFT1(1:n1/2)),'color','blue'); hold on;
xlabel('Frequency f');
ylabel('Fourier transform fft(M)');
title('Fourier transform of Signal')
xlim([0, 25]);
legend('Sampling frequency 500','Sampling frequency 50');

Answer 1

[William Rose]

@Martina Kocan,

The normalization, or vertical scale, of the FFT is affected if you sample the same physical signal at different rates. You can choose to scale the FFT so this does not happen, by diviing the amplitudes of the FFT by the sampling rate. This will work if the signal is sampled for the same duration (i.e. same number of seconds) in both cases.

.

[William Rose]

@Martina Kocan,

I completely agree with you that it would be nice if the FFT ampltiudes were only sensitive to the amplitude of the signal being sampled. But that is not the way the discrete Fourier transform (DFT) is defined. (By the way, it is accurate in this discussion to refer to the DFT, rather than the FFT. The FFT is a clever fast way of computing the DFT. The issue we are discussing affects the DFT, whether it is computed with the FFT algorithm or some other way.)

The amplitude of the DFT (FFT) is proportional to the number of samples. Therefore, if you sample for twice as long at the same sampling frequency, or if you sample for the same duraiton but twice as fast, you will have twice as many data points, and the DFT amplitude will be twice as large. See examples below.

f=6;               %frequency of the signal: same for all records
fs1=40; N1=40;                  %record 1 sampling rate and number of points
t1=(0:N1-1)/fs1;                %time values for record 1
x1=cos(2*pi*f*t1);              %record 1
fs2=80; N2=80;                  %record 2 sampling rate and number of points
t2=(0:N2-1)/fs2;                %time values for record 2
x2=cos(2*pi*f*t2);              %record 2
fs3=40; N3=80;                  %record 3 sampling rate and number of points
t3=(0:N3-1)/fs3;                %time values for record 3
x3=cos(2*pi*f*t3);              %record 3
df1=fs1/N1; fr1=(0:N1-1)*df1;   %frequencies for DFT of x1
df2=fs2/N2; fr2=(0:N2-1)*df2;   %frequencies for DFT of x2
df3=fs3/N3; fr3=(0:N3-1)*df3;   %frequencies for DFT of x3
X1=fft(x1); X2=fft(x2); X3=fft(x3); %compute the DFTs
%next: plot results
figure; subplot(211), plot(t1,x1,'-rx',t2,x2,'-g+',t3,x3,'-bo');
xlabel('Time (s)'); legend('x1','x2','x3');
subplot(212), plot(fr1,abs(X1),'-rx',fr2,abs(X2),'-g+',fr3,abs(X3),'-bo');
xlabel('Frequency (Hz)'); legend('|X1|','|X2|','|X3|');

The plots of the DFT show the full DFT, including the part above the Nyquist frequency, which equals fs/2. The amplitude spectrum above the Nyquist frequency is a flipped copy of the amplitude spectrum below the Nyquist frequency.

Record 2 and record 3 have twice as many points as record 1. Record 2 is sampled twice as fast as record 1, for the same duration. Record three is sampled at the same rate as record 1, for twice the duration.

The DFT amplitudes for records 2 and 3 are twice as big as the DFT amplitude for record 1, even though they all are recordings of the same signal with amplitude 1. This is kind of regrettable but it is what it is. One way to understand it is that it is a consequence of the discrete version of Parseval's theroem.

Answer 2

[Paul]

What is meant by "insensitive?"

I don't see how DFTs could, in general, be the same for different samples of a continuous time domain signal. Recall that the elements of the DFT, which are computed by the fft(), are frequency domain samples of the finite duration sequence. So a related (synonymous?) question would be whether or not the DTFTs of both sequences are the same. If the DTFTs are different, then so will be the DFTs, The DTFT for a finite duration sequence can be computed by freqz().

clear all;
M = 0.1;
c=0.15;
% Signal 1
fs1 = 50; %sampling frequency
rs1 = 1/fs1;  
r1 = 0:rs1:2-rs1; 
L1 = length(r1); %signal length
% Gaussian disribution 
M1 = M * exp(-1/2*(r1/c).^2);

% Signal 2
fs2 = 500; %sampling frequency
rs2 = 1/fs2; 
r2 = 0:rs2:2-rs2; 
L2 = length(r2); %signal length
% Gaussian disribution
M2 = M * exp(-1/2*(r2/c).^2);

[h1,w1] = freqz(M1,1,1024);
[h2,w2] = freqz(M2,1,1024);

figure
plot(w1/rs1/2/pi,abs(h1),w2/rs2/2/pi,abs(h2))
ylabel('abs(DTFT');
xlabel('Frequency (Hz)')

Obviously the blue curve is not the red curve. It doesn't cover the same frequency range. Even with normalzation, it's still not the same, even over its limited frequency range.

figure
plot(w1/rs1/2/pi,abs(h1)/abs(h1(1)/h2(1)),w2/rs2/2/pi,abs(h2))
ylabel('abs(DTFT)');
xlabel('Frequency (Hz)')
xlim([0 50])

We can also view the DFT as a sequence versus index, which is the same as looking at the DTFT versus normalized frequency.

plot(w1/2/pi,abs(h1/abs(h1(1)/h2(2))),w2/2/pi,abs(h2))
ylabel('abs(DTFT)');
xlabel('Normalized Frequency')

Even with the scaling, DFT samples of the red curve won't be samples of the blue curve. To illustrate

figure
plot(w1/2/pi,abs(h1/abs(h1(1)/h2(2))),w2/2/pi,abs(h2))
hold on
hdft1 = fft(M1);
hdft2 = fft(M2);
plot((0:(numel(M2)-1))/numel(M2),abs(hdft2),'o')
plot((0:(numel(M1)-1))/numel(M1),abs(hdft1)/abs(hdft1(1)/hdft2(1)),'o')
ylabel('abs()');
xlabel('Normalized Frequency')
xlim([0 .5])