[Matlab] Jak poprawnie wygenerować widmo sygnału za pomocą funkcji fft?

Question

Od kilku dni mam problem z generacja widma do zadanego sygnału - chciałem wykorzystać matlabowska funkcję fft. Mam przykładowo taki kod : Jak w takim wypadku dokonać analizy do wyznaczenia widma sygnału zadanego - co trzeba dodać lub zmodyfikować - chociaż wolałbym nie modyfikować :). Czy takie...

leon99 · Answer

Sygnał to y

wydaje się poprawne, uruchomiłem to w Matlabie i działa, masz tylko duży przeciek widma (leakage) i trzeba wziąć tylko połowę wykresu, bo jak wiadomo dla sygnałów rzeczywistych , widmo jest powielone

Hunter-Killer · Answer

Tak wiem - tutaj brałem pod uwagę również ujemne wartości - dlatego ten 2 prążek.

Tylko czy widmo nie powinno mieć na osi x częstotliwości ?

Oraz zastanawiam się jak podpisać 2 wykres :

stem(f,Pyy(1:513))

leon99 · Answer

tak, na osi x ma być częstotliwość można nazwać Widmo mocy (power spectrum)

Hunter-Killer · Answer

figure(1)

A= 3; %amplituda sygnału
lo= 3; %liczba okresów
P= 0; %zacznij na czas=0
TE = 6; %zakończ na czas = 6
lp=112; %liczba próbek
o=3; %czas trwania okresu sygnalu
w = pi*lo/TE;
plp=TE/(lp-1); %obliczenie co jaki czas ma pojawić się próbka
m=P:plp:TE;
y=A*sin(2*w*m);
plot(m,y)
Title(['\bf Sygnał zadany'],'Color','black');
grid on
xlabel('Czas T [sek]')
ylabel('Napięcie U [V]')

figure(2)
Y = abs(fft(y,1024));
Pyy = Y.* conj(Y) / 1024;
f = 1000*(0:512)/1024;
stem(f,Pyy(1:513))
Title(['\bf Widmo mocy'],'Color','black');
grid on
xlabel('')
ylabel('')

figure(3)
cze=(1/TE);
g=linspace(P,cze,1024)
stem(g,Y)
Title(['\bf Widmo sygnału zadanego'],'Color','black');
grid on
xlabel('Częstotliwość [Hz]')
ylabel('')

W ten sposób myślę - mam tylko problem z opisaniem osi brakujących.

leon99 · Answer

xlabel na rys 2 i 3 to częstotliwość

nie zgadza mi się skalowanie osi

amplitudę sygnału masz 3 czyli moc powinna być ok. 9 a nie 27 (około bo jest przeciek , czyli brak dopasowania liczby próbek do okresu sinusoidy.

Wzorcowe programy w Matlabie i wiadomości o DFT masz tutaj: (kopiuj z najlepszych :-)

1. Teoria

2. Praktyka

Promuję tematy:
09.11.2009 Ładowanie akumulatorków w mierniku MPI 511 - jak obejść zabezpieczenie?
Mam na sprzedaż:
03.03.2014 Przekaźnik spolaryzowany TGL 6625 10 szt NOWE

Hunter-Killer · Answer

Ja rozumiem tylko, ze ja robię taki program , ze samemu się podaję liczbę próbek - i to jest z rożnym skutkiem - wystarczy, ze dam liczbę 64 i jest tak jak mówisz.

Nie chciałem robić tak by sam program ja dopasowywał.

Hunter-Killer · Answer

Pokazałem te wyniki pewnej wiarygodnej osobie - niby jest praktycznie całe źle - ponieważ, faktycznie maja być 2 prążki ale nie w tych miejscach...

analfabet · Answer

Pobieżnie przeczytałem temat, nie jestem ekspertem w dziedzinie widm, ale nie jest tak, że do szybkiej transformaty Fouriera potrzeba 2^n próbek (i dlatego dla 64 działa)?

Hunter-Killer · Answer

Znaczy ja pokazywałem Widmo sygnału zadanego i tam jest błędnie - a z tymi 64 to było gęstości widma..

mungo · Answer

Wrzucę ci kod jaki na zajęciach z przetwarzania sygnałów używaliśmy do widma, wyciągnięte z większego programu ale powinno działać. Jak coś to pytaj.


Nf=2^11;
N21=Nf/2+1;
w=fft(x,Nf) %x to sygnał
v=abs(w);
Nf=length(v);
tf=linspace(0,Fs/2,N21); %Fs czestotliwosc probkowania
plot(tf,v(1:N21));
xlabel('czestotliwosc[Hz]')
ylabel('modul widma')

Hunter-Killer · Answer

Wydaje się poprawne - tylko zastanawiam się jak to dopasować do tego fragmentu jaki ja mam - bym mógł sobie ładnie modelować sygnał :/

mungo · Answer

Poprawne musi być bo dyktowane przez pewnego dr. A co dokładnie chcesz robić? Co z palca ma być wpisywane do generowania sygnału?

Hunter-Killer · Answer

No ja mam taki większy program - ale idea jest ta sama - podaje się okres sygnału, amplitudę, typ (sinus np), liczbe okresów, od do ile pokazać :

A= 3; %amplituda sygnału
lo= 3; %liczba okresów
P= 0; %zacznij na czas=0
TE = 6; %zakończ na czas = 6
lp=112; %liczba próbek
o=3; %czas trwania okresu sygnalu
w = pi*lo/TE;
plp=TE/(lp-1); %obliczenie co jaki czas ma pojawić się próbka
m=P:plp:TE;
y=A*sin(2*w*m);
plot(m,y)
Title(['\bf Sygnał zadany'],'Color','black');
grid on
xlabel('Czas T [sek]')
ylabel('Napięcie U [V]')

czyli jak powyżej - i do tak zadanego sygnału potrzebuje by wyznaczył mi matlab Widmo sygnału(z dft czy fft)

mungo · Answer

Nie wiem czy dobrze rozumiem, ale jeśli ta twoja część generuje sygnał to dodaj: Fs=1/plp; i dodaj tą moją cześć tylko tam w linijce w=fft(x,Nf); x zamien na y

Hunter-Killer · Answer

Wyszło tak - i mi się wydaję, że jest źle umiejscowiony prążek - zakładając, ze w tej wersji co mi podałeś brany jest pod uwagę 0,5*okres sygnału - bo gdyby był pełny okres to by musiały być 2 prążki

http://imageshack.us/photo/my-images/88/wykresv.jpg

mungo · Answer

Dobrze jest tak powinno być... jeden prążek to znaczy że tylko jest jedna częstotliwość składowa (jeden sinus) i ten prążek jest na wysokości tej częstotliwości...
Spróbuj sobie zrobić np:


Tc=1; %czas próbki
Fs=100; %czestotliwosc próbkowania
t=(0:1/Fs:Tc); %wektor czasu
x=sin(2*pi*15*t)+sin(2*pi*40*t); %t wektor czasu, 15 i 40 częstotliwości sinusów
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%WIDMO

Nf=2^11;
N21=Nf/2+1;
w=fft(x,Nf) %x to sygnał
v=abs(w);
Nf=length(v);
tf=linspace(0,Fs/2,N21); %Fs czestotliwosc probkowania
plot(tf,v(1:N21));
xlabel('czestotliwosc[Hz]')
ylabel('modul widma')

Masz sygnał składający się z sumy 2 sinusów(2 częstotliwości składowe) i powinieneś dostać 2 prążki jeden na 15Hz i 2 na 40Hz

Hunter-Killer · Answer

No właśnie - na wysokości częstotliwości - u mnie 0 - a to na pewno jest źle...

mungo · Answer

Powiększ sobie wykres, musi być dobrze... Przy okazji policz sobie jaką masz częstotliwość tego swojego sygnału bardzo mała będzie.

I tak przy okazji sprawdź co ci pokaże ten mój skrypt z postu wyżej i jak coś zmień Fs na 1000, bo coś koślawy sinus może wyjść. Jak dalej będzie w zerze to jutro zainstaluje sobie matlaba i sprawdzę dokładnie jak to działa bo na razie prawie z głowy pisze.

Hunter-Killer · Answer

Powiększyłem i jest źle - zresztą po skali już widać , ze coś nie pasuje - u Ciebie działa - chyba się nie da przerobić Twojej wersji tak by współgrała z moja - tj. dowolnym modelowaniem wykresu - też się zastanawiam bo z tego co widzę to u Ciebie jest wymagane, za każdym razem dopisywanie próbek 0 - a przecież jeśli bym dokładnie trafił w odpowiednie punkty to nie musiałbym ich dopisywać i bym nie miał przecieku widma - a ja tak spróbowałem, ze idealnie trafia w maximum - a gdy zmniejszyłem liczbę próbek 0 to widmo tez juz nie pasowało.

Hunter-Killer · Answer

Tak się zastanawiam, może lepiej by było wprowadzić (jakoś) z definicji wzór na DFT - jest chyba dokładniejszy i sam bym mógł dokładnie jakoś dopasować go do swoich danych.

Znalazłem DFT - ale coś nie działa - mimo, ze wygląda ok :

function a = mydft(x,N)
%MYDFT Calculates the discrete fourier transform
%usage:
%[a]=mydft(x)
%x=[ x[0] x[1] ... x[N-1] ] - vector containing discrete time sequence
%a=[ a_0 a_1 ... a_N-1 ] - vector of discrete time Fourier Series
%coefficients, k=0,...,N-1
for k=0:N-1
for m=0:N-1
w=exp(-j*2*pi/N.*m.*k);
end
a=sum(x(m).*w);
end
end

mungo · Answer

Dobra zainstalowałem sobie matlaba, można coś pokombinować...

Hunter-Killer napisał:
Wyszło tak - i mi się wydaję, że jest źle umiejscowiony prążek - zakładając, ze w tej wersji co mi podałeś brany jest pod uwagę 0,5*okres sygnału - bo gdyby był pełny okres to by musiały być 2 prążki

http://imageshack.us/photo/my-images/88/wykresv.jpg

Brana jest pod uwagę częstotliwość próbkowania, zgodnie z twierdzeniem Shannona-Kotielnikowa, częstotliwość próbkowania musi być co najmniej 2 razy większa niż częstotliwość sygnału, czyli próbkujemy z częstotliwością np 1000 Hz maksymalna częstotliwość sygnału może wynieść 500Hz i dlatego w tym skrypcie częstotliwość jest wyskalowana do połowy częstotliwości próbkowania czyli 500Hz, jeśli tego by nie było oś biegła by do 1000Hz tylko że od punktu 500Hz było by odbicie lustrzane.

Co do tego zdjęcia wygenerowałem sobie taki sam przebieg i przepuściłem go przez ten mój skrypt do transformaty, wszystko jest w porządku:
próbkowanie 100Hz

Na pierwszy rzut oka wydaje się że jest w zerze.

próbkowanie 10Hz

Skrypt:


Fs=10; % lub Fs=100;
Tc=6;
t=(0:1/Fs:Tc);
x=sin(2*pi*t*1/6);
subplot(211)
plot(t,x);

Nf=2^11;
N21=Nf/2+1;
w=fft(x,Nf);
v=abs(w)
tf=linspace(0,Fs/2,N21);
subplot(212)
plot(tf,v(1:N21))

Zastanowił bym się raczej nad tym twoim skryptem do generowania sygnału bo część od fft jest w pożądku i z tym nic bym nie kombinował.

A i jeszcze coś, podał byś kod jaki użyłeś do tego sygnału ze zdjęcia bo coś mi tam nie pasuje.

Hunter-Killer · Answer

No ja właśnie też sobie przerobiłem program - już nie podaję ilości próbek - tylko częstotliwość próbkowania, oraz nie podaje długości okresu tylko częstotliwość i w ten sposób tworze zadany sygnał i potem widmo i też mam ok - wydaje mi się , ze przy poprzednim zapisie (jak podawałem ilość próbek, długość okresu itd i tak zapisany sygnał przechodził, przez fft w matlabie) nie dało się za bardzo poprawnie zrobić widma sygnału zadanego - nie potrafił pod pasować zmiennych i dlatego nie wychodziło poprawnie (stąd był ten mój screen co podsyłałem) teraz mam tak :

Częstotliwość próbkowania - 16
Amplituda 1 V
Sinusoida 3Hz

http://imageshack.us/photo/my-images/830/wykreszd.jpg/

Jak widać mam najwyższy prążek przy 3Hz - czyli pasuje i jest przeciek, ponieważ też próbkowanie nie było "za gęste".

Elektronik9 · Answer

Przeciek widma można skorygować dobierając odpowiednio liczbę próbek w zadanym przedziale, tak by była ona całkowita. Nie pamiętam dokładnie, zajrzyj do jakiejś książki od DSP, np. do Lyonsa.

Hunter-Killer · Answer

Zgadza się - ja tylko napisałem ,ze jest przeciek i jest to prawidłowe w tym wypadku .

Hunter-Killer · Answer

A to fft z matlaba można wykorzystywać do sygnału piłokształtnych ? Tego typu: x = (A/2)*mod(t*F,T);

[Matlab] Jak poprawnie wygenerować widmo sygnału za pomocą funkcji fft?

Post #1

Post #2

Post #3

Post #4

Post #5

Post #6

Post #7

Post #8

Post #9

Post #10

Post #11

Post #12

Post #13

Post #14

Post #15

Post #16

Post #17

Post #18

Post #19

Post #20

Post #21

Post #22

Post #23

Post #24

Post #25

Post #26

Podsumowanie tematu