Elektroda.pl
Elektroda.pl
X
Proszę, dodaj wyjątek dla www.elektroda.pl do Adblock.
Dzięki temu, że oglądasz reklamy, wspierasz portal i użytkowników.

[Matlab] Widmo do zadanego sygnału

18 Lis 2011 16:25 11418 25
  • Poziom 9  
    Od kilku dni mam problem z generacja widma do zadanego sygnału - chciałem wykorzystać matlabowska funkcję fft. Mam przykładowo taki kod :

    Code:
    figure(1)
    
    A= 4; %amplituda sygnału
    lo= 3; %liczba okresów
    P= 0; %zacznij na czas=0
    TE = 6; %zakończ na czas = 6
    lp=44; %liczba próbek
    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)


    Jak w takim wypadku dokonać analizy do wyznaczenia widma sygnału zadanego - co trzeba dodać lub zmodyfikować - chociaż wolałbym nie modyfikować :).

    Czy takie rozwiązanie jest poprawne ?

    Code:
    figure(2)
    
    Y = abs(fft(y,1024));  %wyznacz wartości z analizy fft
    Pyy = Y.* conj(Y) / 1024;
    f = 1000*(0:512)/1024;
    stem(f,Pyy(1:513))

    figure(3)
    f=linspace(P,TE,1024)
    stem(f,Y)
  • Poziom 19  
    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
  • Poziom 9  
    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))
  • Poziom 9  
    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.
  • Poziom 19  
    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
  • Poziom 9  
    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ł.
  • Poziom 9  
    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...
  • Poziom 12  
    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)?
  • Poziom 9  
    Znaczy ja pokazywałem "Widmo sygnału zadanego" i tam jest błędnie - a z tymi 64 to było gęstości widma..
  • Poziom 19  
    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.
    Code:

    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')
  • Poziom 9  
    Wydaje się poprawne - tylko zastanawiam się jak to dopasować do tego fragmentu jaki ja mam - bym mógł sobie ładnie modelować sygnał :/
  • Poziom 19  
    Poprawne musi być bo dyktowane przez pewnego dr.
    A co dokładnie chcesz robić? Co z palca ma być wpisywane do generowania sygnału?
  • Poziom 9  
    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)
  • Poziom 19  
    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
  • Poziom 9  
    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
  • Poziom 19  
    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:
    Code:

    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
  • Poziom 9  
    No właśnie - na wysokości częstotliwości - u mnie 0 - a to na pewno jest źle...
  • Poziom 19  
    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.
  • Poziom 9  
    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.
  • Poziom 9  
    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
  • Poziom 19  
    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
    [Matlab] Widmo do zadanego sygnału
    Na pierwszy rzut oka wydaje się że jest w zerze.

    próbkowanie 10Hz
    [Matlab] Widmo do zadanego sygnału
    Skrypt:
    Code:

    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.
  • Poziom 9  
    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".
  • Poziom 30  
    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.
  • Poziom 9  
    Zgadza się - ja tylko napisałem ,ze jest przeciek i jest to prawidłowe w tym wypadku .
  • Poziom 9  
    A to fft z matlaba można wykorzystywać do sygnału piłokształtnych ?

    Tego typu:

    x = (A/2)*mod(t*F,T);