FFT of NMR signal: frequency cut off

11 visualizzazioni (ultimi 30 giorni)
Nicolai
Nicolai il 28 Giu 2014
Risposto: Aaron Globisch il 15 Mag 2017
Hello,
I've got a problem with the fourier transform of a time domain nmr signal. When I perform the fourier transform I receive a weired spectrum. The first Peak is cutted in half, as you can see in the attached picture. Does somebody know, why lower frequencies are not preserved? Is there a mistake in my code or does somebody have an idea? The raw data is ok and contain all informations as I know when I process the data with standard nmr software.
Many thanks, N.
clear all;
close all;
clc
tdomain = load('nmr_ethanol.mat');
tdomain = tdomain.signal;
AqT = tdomain(1,end);
SampRate = tdomain(2,1)-tdomain(1,1);
maxfreq = (length(tdomain(:,1))-1)./(length(tdomain(:,1))*SampRate);
lowfrq = -2295.07256526605;
FreqVek = linspace(0,maxfreq,length(tdomain(:,1)))+(lowfrq);
FreqVek2 = linspace(0,maxfreq,2*length(tdomain(:,1)))+(lowfrq);
npoints = 2.*length(tdomain(:,1));
% fourier transformation
fdomain = fft(tdomain(:,2),npoints);
% cut off symmetrical part and frequency shift
fdomain_sym = fdomain(1:npoints./2);
fdomain_shift = fftshift(fdomain_sym/npoints);
% Plot
fig1 = figure(1);
subplot(2,2,1)
plot(tdomain(:,1),tdomain(:,2));
title('time domain')
xlabel('time, t / s');
xlim([0 0.75])
subplot(2,2,2)
plot(FreqVek2,real(fdomain),'-');
grid on;
title('frequency domain')
xlabel('frequency, f / s^{-1}');
xlim([-2500 3300]);
subplot(2,2,3)
plot(FreqVek,real(fdomain_sym),'-');
grid on;
title('frequency domain, symmetrical half')
xlabel('frequency, f / s^{-1}');
xlim([-2500 -1500]);
subplot(2,2,4)
plot(FreqVek,real(fdomain_shift),'-');
grid on;
title('frequency domain, frequency shift')
xlabel('frequency, f / s^{-1}');
xlim([200 1000]);

Risposte (4)

Nicolai
Nicolai il 29 Giu 2014
Thanks a lot for your answers. Unfortunatly the change in the plotting settings do not solve the problem. In nmr spectroscopy we always look at the real part of the spectrum and do not use the magnitude spectrum. But thanks for the suggestions!
Could there be a reason why I cannot restore all frequencies?
The spectrum supposed to be look like this:
Thanks a lot for all your ideas! N.

Star Strider
Star Strider il 28 Giu 2014
Not certain about this because I don’t have your data, but see if changing your plot statements in the last two supblot calls from real to abs:
plot(FreqVek,abs(fdomain_sym),'-');
and
plot(FreqVek,abs(fdomain_shift),'-');
solves your problem.

Daniel kiracofe
Daniel kiracofe il 28 Giu 2014
Hi Nicolai. I'm not sure exactly what you are expecting to see, but I suspect you just didn't multiply by 2 when converting from 2 sided to 1 sided spectrum. Generally, I never use matlab's fft() directly. I always use a wrapper script. Below is a basic wrapper script. You can find more discussion of it at my website <http://www.mechanicalvibration.com/Making_matlab_s_fft_functio.html>
% function [frq, amp, phase] = simpleFFT( signal, ScanRate)
% Purpose: perform an FFT of a real-valued input signal, and generate the single-sided
% output, in amplitude and phase, scaled to the same units as the input.
%inputs:
% signal: the signal to transform
% ScanRate: the sampling frequency (in Hertz)
% outputs:
% frq: a vector of frequency points (in Hertz)
% amp: a vector of amplitudes (same units as the input signal)
% phase: a vector of phases (in radians)
function [frq, amp, phase] = simpleFFT( signal, ScanRate)
n = length(signal);
z = fft(signal, n); %do the actual work
%generate the vector of frequencies
halfn = n / 2;
deltaf = 1 / ( n / ScanRate);
frq = (0:(halfn-1)) * deltaf;
% convert from 2 sided spectrum to 1 sided
%(assuming that the input is a real signal)
amp(1) = abs(z(1)) ./ (n);
amp(2:halfn) = abs(z(2:halfn)) ./ (n / 2);
phase = angle(z(1:halfn));

Aaron Globisch
Aaron Globisch il 15 Mag 2017
Hello,
i figured out the correct way:
specturm = abs(fftshift(fft(data, s_points)/s_points));
plot(FreqVek2, spectrum);
where s_points is the amount of sampled points and FreqVek2 the frequency vector just like u declared it.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by