How to resample frequency domain signal?

Hi!

I have a spectrum with freqency resolution 0.01875 Hz, and I want to have 1Hz. Is there a way to convert the signal to have 1Hz frequency resolution without deforming and ruining either the phase or the amplitude of the signal? (I can not change the length or the sampling frequency of the time domain signal).

[Mathieu NOE]

hello

if you have magnitude and phase as vectors (vs frequency) , you can simply use interp1 with a new frequency vector ( delta = 1) to get the same graph with a different frequency resolution

[Csanad Levente Balogh]

Thank you! But does this account for the energy of the spectrum. I think, by getting rid of some frequncies, you lower the energy stored in the signal. Is there a way to keep it the same?

[Mathieu NOE]

ok - so you would a kind of energy (rms) averaging over the original samples within each 1 Hz band ?

are you willing to share your data ?

[Csanad Levente Balogh]

In data you can find a vector containing the frequencies and another vector containing the spectrum values. Single sided spectrum is calculated like:

x time domain signal
L length of time domain signal
fs sampling frequency

Y = fft(x);
Y = Y/L;
Y = Y(1:L/2+1);
Y(2:end-1) = 2*Y(2:end-1);

everywhere in the code, and the frequency is calculated like:

f = (fs/L)*(0:(L/2));

My problem is, I want to reduce the number of frequency components in the spectrum, without ruinning the amplitude of phase, and without lowering the energy stored in the signal. (altough I know that I necessarily loose information by reducing the number of frequeny components).

Answer 1

[Mathieu NOE]

hello

• first code below : this is a script showing true averaging (even on complex Y data) over N samples

your actual frequency resolution is 0.0625 Hz, so to get results with 1 Hz spacing you need to do N = 16 samples averaging

as Y (the fft data ) is complex , you can see that summing complex data can lead to "apparent" loss of magnitude, this is because the phase changes very abruptly from one frequency to the next

so summation of 2 complex with almost opposite phase can lead to much lower magnitude in the summation;

but this preserves the best the phase information

• the second code below acts like a "autospectrum" averaging , we average the magnitude and the phase separately, so it better preserve the magnitude , but the resulting angle is not really meaningfull to me ; by the way "autospectrum" averaging is never displayed for phase - only magnitude

Code 1

load('data.mat')

N = 16;             % Number of elements to create the mean over

%% slow for loop method 
% zero overlap mean averaging
[m,n] = size(f);

for ci=1:fix(m/ N)
    start_index = 1+(ci-1)*N;
    stop_index = min(start_index+ N -1,m);
    avg_f(ci) =mean(f(start_index:stop_index));  % 
    avg_Y(ci) =mean(Y(start_index:stop_index));  % 
end

%% faster method 
f_new = my_mean_N_samples(f,N);
Y_new = my_mean_N_samples(Y,N);


%% plot
subplot(2,1,1),semilogy(f, abs(Y),'b',f_new, abs(Y_new),'dr',avg_f, abs(avg_Y),'g');
legend('original data','averaged with slow for loop','averaged with faster method');
xlim([3950 4050]);
xlabel('Hz');
ylabel('Magnitude');

subplot(2,1,2),plot(f, angle(Y),'b',f_new, angle(Y_new),'dr',avg_f, angle(avg_Y),'g');
legend('original data','averaged with slow for loop','averaged with faster method');
xlim([3950 4050]);
xlabel('Hz');
ylabel('Phase');

function out = my_mean_N_samples(x,n)
    x = x(:);             % make sure the data are column oriented
    s1 = size(x, 1);      % Find the next smaller multiple of n
    m  = s1 - mod(s1, n);
    y  = reshape(x(1:m), n, []);     % Reshape x to a [n, m/n] matrix
    out = transpose(sum(y, 1) / n);  % Calculate the mean over the 1st dimension

end

code 2 :

load('data.mat')


N = 16;             % Number of elements to create the mean over

Ymod = abs(Y);
Yangl = angle(Y);


%% slow for loop method 
% zero overlap mean averaging
[m,n] = size(f);

for ci=1:fix(m/ N)
    start_index = 1+(ci-1)*N;
    stop_index = min(start_index+ N -1,m);
    avg_f(ci) =mean(f(start_index:stop_index));  % 
    avg_Ymod(ci) =mean(Ymod(start_index:stop_index));  % 
    avg_Yangl(ci) =mean(Yangl(start_index:stop_index));  % 
end

%% faster method 
f_new = my_mean_N_samples(f,N);
Ymod_new = my_mean_N_samples(Ymod,N);
Yangl_new = my_mean_N_samples(Yangl,N);



%% plot
subplot(2,1,1),semilogy(f, Ymod,'b',f_new, Ymod_new,'dr',avg_f, avg_Ymod,'g');
legend('original data','averaged with slow for loop','averaged with faster method');
xlim([3950 4050]);
xlabel('Hz');
ylabel('Magnitude');

subplot(2,1,2),plot(f, Yangl,'b',f_new, Yangl_new,'dr',avg_f, avg_Yangl,'g');
legend('original data','averaged with slow for loop','averaged with faster method');
xlim([3950 4050]);
xlabel('Hz');
ylabel('Phase');



function out = my_mean_N_samples(x,n)
    x = x(:);             % make sure the data are column oriented
    s1 = size(x, 1);      % Find the next smaller multiple of n
    m  = s1 - mod(s1, n);
    y  = reshape(x(1:m), n, []);     % Reshape x to a [n, m/n] matrix
    out = transpose(sum(y, 1) / n);  % Calculate the mean over the 1st dimension

end