MATLAB FFT OśX granice brudząc i fftshift

Question

jest to pierwszy raz, kiedy m użyciu funkcji FFT i Próbuję wykreślić widma częstotliwości prostej funkcji cosinus:

F = cos (2 * pi * 300 * t)

Częstotliwość próbkowania to 220500. Wykreślam jedną sekundę funkcji f.

Oto moja próba:

time = 1; 
freq = 220500; 
t = 0 : 1/freq : 1 - 1/freq; 
N = length(t); 
df = freq/(N*time); 

F = fftshift(fft(cos(2*pi*300*t))/N); 
faxis = -N/2/time : df : (N/2-1)/time; 

plot(faxis, real(F)); 
grid on; 
xlim([-500, 500]); 

Dlaczego mam dziwne wyniki, kiedy zwiększyć częstotliwość do 900Hz? Te nieparzyste wyniki można naprawić, zwiększając ograniczenia osi X z, powiedzmy, od 500 Hz do 1000 Hz. Czy to jest prawidłowe podejście? Zauważyłem, że wielu innych ludzi nie używało fftshift(X) (ale myślę, że oni tylko zrobili jednostronną analizę widma).

Dziękuję.

Answer 1

Oto moja odpowiedź zgodnie z obietnicą.

Pierwsze lub twoje pytania dotyczące tego, dlaczego "uzyskujesz dziwne wyniki, gdy zwiększasz częstotliwość do 900 Hz" są powiązane z funkcją skalowania wykresu Matlab opisaną przez @Castilho. Gdy zmienisz zakres osi X, Matlab będzie starał się być pomocny i przeskalować oś y. Jeśli wartości szczytowe wykraczają poza określony zakres, matlab przybliży małe błędy numeryczne generowane w procesie. Możesz to naprawić poleceniem "ylim", jeśli ci to przeszkadza.

Jednak twoje drugie, bardziej otwarte pytanie "czy to właściwe podejście?" wymaga głębszej dyskusji. Pozwólcie, że powiem wam, w jaki sposób chciałbym zastosować bardziej elastyczne rozwiązanie, aby osiągnąć swój cel polegający na wykreślaniu fali cosinus.

Zaczynasz z następujących czynności:

time = 1; 
freq = 220500; 

Ten natychmiast podnosi alarm w mojej głowie. Patrząc na resztę postu, wydaje się, że jesteś zainteresowany częstotliwościami w zakresie sub-kHz. Jeśli tak jest, to częstotliwość próbkowania jest zbyt duża, ponieważ limit Nyquista (sr/2) dla tej szybkości przekracza 100 kHz.Zgaduję, że zamierzałeś użyć wspólnej częstotliwości próbkowania audio 22050 Hz (ale mogę się mylić tutaj)?

Tak czy inaczej, twoja analiza kończy się numerycznie OK. Jednak nie pomaga sobie zrozumieć, w jaki sposób FFT można najskuteczniej wykorzystać do analizy w rzeczywistych sytuacjach.

Pozwolę sobie opublikować, jak to zrobię. Poniższy skrypt wykonuje prawie dokładnie to, co robi twój skrypt, ale otwiera pewien potencjał, na którym możemy go zbudować. .

%// These are the user parameters 
durT = 1; 
fs = 22050; 
NFFT = durT*fs; 
sigFreq = 300; 

%//Calculate time axis 
dt = 1/fs; 
tAxis = 0:dt:(durT-dt); 

%//Calculate frequency axis 
df = fs/NFFT; 
fAxis = 0:df:(fs-df); 

%//Calculate time domain signal and convert to frequency domain 
x = cos( 2*pi*sigFreq*tAxis ); 
F = abs( fft(x, NFFT)/NFFT ); 

subplot(2,1,1); 
plot( fAxis, 2*F ) 
xlim([0 2*sigFreq]) 
title('single sided spectrum') 

subplot(2,1,2); 
plot( fAxis-fs/2, fftshift(F) ) 
xlim([-2*sigFreq 2*sigFreq]) 
title('whole fft-shifted spectrum') 

Użytkownik oblicza oś czasu i oblicza liczbę punktów FFT z długości osi czasu. To bardzo dziwne. Problem z tym podejściem polega na tym, że rozdzielczość częstotliwościowa fft zmienia się wraz ze zmianą czasu trwania sygnału wejściowego, ponieważ N zależy od zmiennej "time". Polecenie matlab fft użyje rozmiaru FFT odpowiadającego rozmiarowi sygnału wejściowego.

W moim przykładzie obliczam oś częstotliwości bezpośrednio z NFFT. Jest to nieco nieistotne w kontekście powyższego przykładu, ponieważ ustawiam NFFT tak, aby był równy liczbie próbek w sygnale. Jednak użycie tego formatu pomaga zdemistyfikować twoje myślenie i staje się bardzo ważne w moim następnym przykładzie.

** BOCZNA UWAGA: W swoim przykładzie używasz rzeczywistego (F). O ile nie masz bardzo dobrego powodu, aby wyodrębnić tylko prawdziwą część wyniku FFT, znacznie częściej można wyodrębnić wielkość FFT używając abs (F). Jest to odpowiednik sqrt (rzeczywisty (F).^2 + imag (F).^2). **

W większości przypadków będziesz chciał użyć krótszego NFFT. Może to być spowodowane tym, że prawdopodobnie przeprowadzasz analizę w systemie czasu rzeczywistego, lub ponieważ chcesz uśrednić wynik wielu FFT razem, aby uzyskać pojęcie o spektrum średniej dla sygnału zmieniającego się w czasie, lub ponieważ chcesz porównać widma sygnały o różnym czasie trwania bez marnowania informacji. Po prostu za pomocą polecenia fft o wartości NFFT < liczba elementów w twoim sygnale spowoduje, że fft zostanie obliczony na podstawie ostatnich punktów NFFT sygnału. To trochę marnotrawstwo.

Poniższy przykład jest bardziej odpowiedni dla użytecznej aplikacji. To pokazuje, w jaki sposób rozdzielić sygnał na bloki, a następnie przetworzyć każdy blok i uśrednić wynik:

%//These are the user parameters 
durT = 1; 
fs = 22050; 
NFFT = 2048; 
sigFreq = 300; 

%//Calculate time axis 
dt = 1/fs; 
tAxis = dt:dt:(durT-dt); 

%//Calculate frequency axis 
df = fs/NFFT; 
fAxis = 0:df:(fs-df); 

%//Calculate time domain signal 
x = cos( 2*pi*sigFreq*tAxis ); 

%//Buffer it and window 
win = hamming(NFFT);%//chose window type based on your application 
x = buffer(x, NFFT, NFFT/2); %// 50% overlap between frames in this instance 
x = x(:, 2:end-1); %//optional step to remove zero padded frames 
x = ( x' * diag(win) )'; %//efficiently window each frame using matrix algebra 

%// Calculate mean FFT 
F = abs( fft(x, NFFT)/sum(win) ); 
F = mean(F,2); 

subplot(2,1,1); 
plot( fAxis, 2*F ) 
xlim([0 2*sigFreq]) 
title('single sided spectrum') 

subplot(2,1,2); 
plot( fAxis-fs/2, fftshift(F) ) 
xlim([-2*sigFreq 2*sigFreq]) 
title('whole fft-shifted spectrum') 

używam okno Hamminga w powyższym przykładzie. Okno, które wybierzesz, powinno pasować do wybranej aplikacji. Wybrana wielkość zależy w pewnym stopniu od typu używanego okna. W powyższym przykładzie okno Hamminga obciąża próbki w każdym buforze w kierunku zera od środka każdej ramki. Aby wykorzystać wszystkie informacje w sygnale wejściowym, ważne jest, aby niektóre nakładały się. Jeśli jednak używasz prostego prostokątnego okna, nakładanie staje się bezcelowe, ponieważ wszystkie próbki są ważone równo. Im więcej nakładania się używasz, tym więcej przetwarzania jest wymagane do obliczenia średniego widma.

Mam nadzieję, że to pomaga zrozumieć.

Answer 2

Twój wynik jest całkowicie poprawny. Również obliczanie osi częstotliwości jest prawidłowe. Problem leży na skali osi y. Kiedy używasz funkcji xlims, matlab automatycznie przelicza skalę y, dzięki czemu możesz zobaczyć "znaczące" dane. Kiedy pik cosinusa znajduje się poza wyznaczonym przez ciebie limitem (przy f> 500Hz), nie ma żadnych pików do pokazania, więc skala jest obliczana na podstawie niewielkiego szumu (tutaj przy moim komputerze, z matlab 2011a, skala y wynosiła 10 -16).

Zmiana limitu jest rzeczywiście poprawnym podejściem, ponieważ jeśli go nie zmienisz, nie zobaczysz pików na spektrum częstotliwości.

Jedno jednak zauważyłem. Czy jest jakiś powód, abyś wykreślił prawdziwą część transformacji? Zwykle jest to abs(F), które zostanie naniesione na wykres, a nie na rzeczywistą część.

edytuj: W rzeczywistości twoja osi częstotliwości jest słuszna tylko dlatego, że df, w tym przypadku, wynosi 1. Wiersz faxis ma rację, ale obliczenia df nie są prawidłowe.

FFT oblicza N punktów od -Fs/2 do Fs/2. Zatem N punktów w zakresie Fs daje df Fs/N. Jako N/czas = Fs => czas = N/Fs. Zastępując to wyrażeniem df, którego użyłeś: your_df = Fs/N * (N/Fs) = (Fs/N)^2. Jako Fs/N = 1, ostateczny wynik był prawidłowy: P