kwiecień 28th, 2009 marcin
Autorem tego wpisu jest Łukasz Ambroziński, a sam wpis uzyskał pierwsze miejsce w konkursie na najlepszy wpis.
Wielu początkujących gitarzystów boryka się z problemem nastrojenia instrumentu. Jakiś czas temu najlepszą metodą było zwrócenie się do bardziej doświadczonego kolegi, dziś z pomocą mogą przyjść bardzo popularne stroiki elektroniczne, ale także Matlab…
W celu stworzenia odpowiedniego narzędzia będzie nam potrzebna pewna wiedza z teorii muzyki – kilka bardzo podstawowych informacji, aby je zdobyć sięgamy do Wikipedii:
Dźwięk - wrażenie słuchowe spowodowane falą akustyczną rozchodzącą się w ośrodku sprężystym.
Jeśli chodzi o parametry tej fali to nas będzie interesować jej częstotliwość – odpowiada ona wysokości dźwięku. Tworzony stroik będzie działał wykorzystując pomiar częstotliwości zarejestrowanego dźwięku. Nastrojenie instrumentu wiąże się z takim dostosowaniem naciągu struny, aby drgała ona z odpowiednią częstotliwością.
Kolejna definicja niezbędna do zrozumienia działania stroika to: półton - czyli najmniejsza odległość między dwoma dźwiękami w systemie temperowanym. Przykładowo dźwięk pustej struny gitary oraz przyciśniętej na pierwszym progu oddalony jest o pół tonu. Dla nas ważna jest informacja, że iloraz częstotliwości tych dźwięków wynosi
np.:

gdzie:
f1 - częstotliwość drgania struny przyciśniętej na pierwszym progu,
f0 - częstotliwość drgania struny pustej.
Można napisać, że:

analogicznie:

czyli:

Powyższe równania prowadzą do wniosku, że częstotliwości kolejnych dźwięków tworzą ciąg geometryczny o ilorazie
, więc aby znaleźć częstotliwość konkretnego dźwięku wystarczy znać częstotliwość dźwięku podstawowego oraz odległość pomiędzy nimi.
Są to w zasadzie wszystkie informacje potrzebne do stworzenia stroika. W celu lepszego poznania natury dźwięku od strony fizyki polecam wszystkim zainteresowanym stronę http://www.fizykon.org/akustyka/akustyka.htm.
Możemy rozpocząć pisanie kodu. W pierwszym kroku obliczymy z jaką częstotliwością powinny drgać struny. Następnie będziemy z obliczoną częstotliwością wzorcową porównywać częstotliwość zmierzoną.
-
a_podst=440; %[Hz] częstotliwość referencyjna
Jest to częstotliwość, do której zgodnie z przyjętymi umowami zazwyczaj stroimy instrumenty. Z taką częstotliwością powinna drgać najcieńsza struna gitary (e) przyciśnięta na 5-tym progu. Od tej częstotliwości obliczymy wartości wymaganych częstotliwości dla kolejnych pustych strun.
Iloraz częstotliwości odpowiadający półtonowi. Jak wspomniałem wcześniej częstotliwości dźwięku tworzą ciąg geometryczny o ilorazie
. Zapamiętujemy tę wartość w zmiennej delta.
-
struny=[0 0 0 0 0 0];
-
struny(1)=a_podst/delta^5; %e
-
struny(2)=a_podst/delta^10;%h
-
struny(3)=a_podst/delta^14;%g
-
struny(4)=a_podst/delta^19;%d
-
struny(5)=a_podst/delta^24;%a
-
struny(6)=a_podst/delta^29;%e1
-
nazwy_strun=’ehgdae‘;
Powyższy fragment kodu tworzy tablicę z wartościami pożądanymi częstotliwości dla kolejnych strun – obliczamy konkretne, interesujące nas wartości elementów ciągu, o którym wspomniałem wcześniej. Znajdowane są one na podstawie częstotliwości referencyjnej oraz odległości w półtonach od dźwięku a.
e – 5 półtonów (progów) niżej niż a,
h – 5 półtonów niżej niż e czyli 10 półtonów niżej niż a itd.
Bardziej zaawansowani muzycy mogą zmieniając wykładniki potęg w powyższym kodzie zmieniać tonację, w której stroją instrument – aby nastroić jedną tonację niżej należy zwiększyć wykładniki o 1, aby nastroić wyżej należy wykładniki zmniejszyć. Można również zastosować strój inny niż klasyczny.
Przykładowe wartości wykładników dla stroju open E: 5, 10 , 13, 17, 22, 29.
Przejdźmy do kolejnego kroku czyli pomiar i analiza sygnału dźwiękowego:
-
ai = analoginput(‘winsound‘);
-
addchannel(ai,1:2);
-
rate=8192; time=1;
-
set(ai,’
SampleRate‘,rate,’
SamplesPerTrigger‘,rate*time
)
Wykorzystujemy tutaj opisany wcześniej na blogu spób użycia karty dźwiękowej jako karty akwizycji danych. Stosujemy rejestrację z częstotliwością próbkowania 8192Hz, czas rejestracji 1s.
-
start(ai);
-
y= getdata(ai);
-
N = 8192;
-
Rejestrujemy fragment dźwięku, a następnie na otrzymanym sygnale wykonujemy szybką transformatę Fouriera, która pozwoli nam otrzymać widmo częstotliwościowe nagranego dźwięku. 8192 punktowa transformata pozwoli nam dla sygnału próbkowanego 8192Hz otrzymać rozdzielczość 1Hz.
-
-
f = rate/N*(0:N-1);
-
f=f(10:N/10);
-
X=X(10:N/10);
-
-
Szukamy częstotliwości odpowiadającej największemu prążkowi na widmie – jest to podstawowa częstotliwość, z jaką drga struna.
-
roznica=10000;
-
-
roztmp=
abs(fszukane - struny
(i));
-
if roztmp
-
roznica=roztmp;
-
-
end
-
end
Powyższa pętla porównuje otrzymany wynik pomiaru z pożądanymi częstotliwościami.
Ponieważ nie wskazujemy programowi która struna jest aktualnie strojona szuka on w tablicy częstotliwości wzorcowych wartości najbardziej zbliżonej do pomiaru, indeks pod którym występuje ta wartość zapisywany jest w zmiennej dźwięk. Jeśli wartość indeksu zgadza się z numerem dostrajanej struny oznacza to, że zbliżamy się z jej naciągiem do odpowiedniej wartości.
-
roznica = fszukane - struny(dzwiek);
-
-
if roznica > 1
-
napis=’<‘;
-
if roznica > 10
-
napis=’<<‘;
-
end
-
elseif roznica<-1
-
napis=’>‘;
-
if roznica<-10
-
napis=’>>‘;
-
end
-
else
-
napis=’OK‘
-
end
W kodzie powyżej generowane są komunikaty dla użytkownika pokazujące czy należy zwiększyć czy zmniejszyć naciąg struny.
W ten sposób powstała całkiem funkcjonalna aplikacja. Wystarczy dorobić GUI i można się cieszyć własnoręcznie zaprogramowanym stroikiem gitarowym. Może on wyglądać następująco:

Posted in MATLAB | 2 komentarzy »
lipiec 9th, 2008 jaromir
Witam.
Tym razem, zgodnie z zapowiedzią, przedstawię akwizycję obrazów z kamer analogowych w sposób programowy. Przypominam, że obsługę kamer zapewnia moduł Image Acquisition Toolbox.
W przykładzie wykorzystana została popularna kamera internetowa USB.
Do zidentyfikowania urządzeń dostępnych w systemie służy polecenie IMAQHWINFO:
-
info = imaqhwinfo
-
info.InstalledAdaptors
Dla kamer USB odpowiedni sterownik to “winvideo“.
W celu skonfigurowania akwizycji potrzebne są jeszcze informacje o obsługiwanych trybach przez dane urządzenie (m.in. rozdzielczość). Znając nazwę urządzenia można to zrealizować przy pomocy tego samego polecenia co poprzednio:
-
info = imaqhwinfo(‘winvideo‘)
-
info.DeviceInfo
-
info.DeviceInfo.SupportedFormats
Aby dokonać akwizycji konieczna jest konfiguracja sposobu pobierania obrazu. W pierwszej kolejności tworzony jest obiekt akwizycji:
-
OBJ = VIDEOINPUT
(ADAPTORNAME,DEVICEID,
FORMAT)
Gdzie: ADAPTORNAME - nazwa urządzenia, DEVICEID - numer urządzenia (karta akwizycji może umożliwiać podłączenie kilku kamer), FORMAT - jeden z obsługiwanych formatów.
Przykładowo:
-
vidobj = videoinput(‘winvideo‘,1,’RGB24_640×480‘)
Mając utworzony obiekt akwizycji możliwa jest inspekcja jego własności lub podgląd obrazu:
-
inspect(vidobj) % własności obiektu akwizycji
-
preview(vidobj) % podgląd obrazu
-
-
closepreview(vidobj) %zamknięcie podglądu
Konfiguracji obiektu akwizycji dokonuje się następującymi poleceniami:
-
triggerconfig(vidobj,’manual‘) % sposób wyzwalania (tutaj ręczny czyli wyzwalany odpowiednim poleceniem)
-
set(vidobj,’
FramesPerTrigger‘,
1) % ilość ramek pobieranych jednorazowo
-
start(vidobj) % uruchomienie obiektu akwizycji
Akwizycja obrazu:
-
trigger(vidobj); % uruchomienie akwizycji
-
pause(.
5);
% opóźnienie programu aby akwizycja została zrealizowana
-
im = getdata(vidobj,1); % pobranie danych z obiektu akwizycji
W powyższym przykładzie wprowadzone zostało niewielkie opóźnienie ze względu na to, iż uruchomienie akwizycji zajmuje pewną ilość czasu. Alternatywnym sposobem jest sprawdzanie (np. w pętli) czy ramka obrazu została pobrana do bufora. Można to sprawdzić poleceniem:
Po zakończeniu akwizycji zalecane jest usunięcie utworzonych obiektów i zwolnienie urządzenia:
Pobrany obraz można zwizualizować poleceniem:
Ewentualnie, do wizualizacji, można wykorzystać modularne narzędzia graficzne dostarczane przez moduł Image Acquisition Toolbox:

Narzędzia te umożliwiają m.in.:
- podgląd obrazu
- podgląd pikseli obrazu w danym regionie (Tools>Pixel Region)
- pomiar odległości na obrazie(Tools>Measure Distance)
- “wycięcie” fragmentu obrazu(Tools>Crop Image)
- poprawę kontrastu w przypadku obrazów w odcieniach szarości lub indeksowanych (Tools>Adjust Contrast)
Posted in MATLAB | 3 komentarzy »
lipiec 1st, 2008 jaromir
Ostatnio na blogu pojawił się wpis pokazujący w jaki sposób dokonać akwizycji sygnału z karty dźwiękowej przy pomocy Data Acquisition Toolbox. W środowisku MATLAB/Simulink możliwa jest również akwizycja obrazów z kamer analogowych oraz cyfrowych. Do tego celu służy Image Acquisition Toolbox.
Z jakimi kamerami potrafi współpracować MATLAB? Dla systemu Windows, w zasadzie wszystkie kamery które wspierają sterowniki WDM (ang. Windows Driver Model) są obsługiwane (m.in. popularne kamery internetowe). Pełną listę obsługiwanych kamer można znaleźć na stronach producenta (www.mathworks.com).
Aby pobrać obrazy z kamery do MATLABa, należy w pierwszej kolejności podłączyć kamerę oraz zainstalować wymagane sterowniki. Ponieważ etap ten zależy od oprogramowania dostarczanego z kamerę, nie będę go opisywał. Tym niemniej, przed przystąpieniem do pracy w MATLABie dobrze jest sprawdzić czy instalacja sterowników kamery przebiegła poprawnie. Można to zrobić albo używając oprogramowania dostarczanego z kamerą, albo pobrać OpenSource’owe oprogramowanie takie jak: AmCap lub VirtualDub. Z reguły popularne kamery USB po podłączeniu do komputera bez problemu pozwalają na podgląd i akwizycje obrazu. Może jednak się zdarzyć, że dysponujemy kamerą analogową posiadającą wyjście S-Video lub Composite. W takim przypadku konieczne jest jeszcze urządzenie nazywane framegrabberem, dokonujące konwersji sygnału analogowego z kamery na sygnał cyfrowy wprowadzany do komputera. Dla osób, które po raz pierwszy mają do czynienia z kamerami i przetwarzaniem obrazów, wyjaśniam że urządzenie to nie jest bardzo tajemnicze i trudno dostępne. Każda karta telewizyjna, pozwalająca na oglądanie TV w komputerze, jest w zasadzie takim framegrabberem. Ponieważ karty takie posiadają dwa wejścia (s-video i composite) oraz tuner TV, mamy trzy źródła sygnału video. Aby móc pobierać obraz z kamery należy ustawić odpowiednie źródło sygnału, odpowiadające podłączeniu kamery. Oprogramowanie dostarczane z kartą TV (oraz AmCap i VirtualDub) posiada opcje VideoCrossbar, pozwalającą na wybór źródła sygnału.
Zakładając, że podłączenie i sprawdzenie kamery powiodło się, przejdę do omówienia akwizycji obrazu w MATLABie. Opis będzie dotyczył wersji R2008a, w której dostępny jest graficzny interfejs pozwalający na konfigurację urządzenia akwizycji obrazów. Programowa akwizycja również jest możliwa (opis zamieszczę następnym razem).
Uruchamiamy MATLABa i w linii poleceń wywołujemy narzędzie IMAQTOOL:
Narzędzie to pozwala w wygodny sposób na konfigurowanie dostępnych w systemie urządzeń akwizycji sygnału video oraz na podgląd obrazu i jego akwizycję do środowiska MATLAB.

Po lewej stronie dostępna jest przeglądarka dostępnych w systemie urządzeń oraz możliwych trybów pracy (Hardware Browser). Po przeciwnej stronie znajduje się okno (Desktop Help) zwierająca pomoc do aktualnie wybranych opcji. Na środku jest okno podglądu obrazu (Preview) oraz poniżej niego zakładki zawierające opcje konfiguracyjne wybranego urządzenia (Acquisition Parameters). W lewym dolnym rogu znajdują się informacje podsumowujące ustawioną konfigurację urządzenia (Information).
Aby skonfigurować urządzenie do akwizycji obrazu należy w oknie Hardware Browser wybrać urządzenie oraz jeden z udostępnianych trybów dla tego urządzenia (tryby mogą być różne dla różnych kamer i sterowników). Tryb pracy definiuje rozdzielczość akwizycji oraz wykorzystywaną przestrzeń barw (np. RGB lub YUV). Informacje na temat przestrzeni barw można znaleźć w internecie - na potrzeby tego blogu wystarczający będzie wybór przestrzeni RGB24, w której każda składowa koloru jest zapisana na 8 bitach (co daje 256 możliwych wartości dla jednej składowej).
Po wybraniu trybu dostępne stają się parametry w oknie Acquisition Parameters. Nie będę ich tutaj omawiał szczegółowo - ich opis jest w dokumentacji (wyświetlanej z prawej strony). Omówię tylko parametry potrzebne do pobrania jednej ramki obrazu.
Możliwa jest akwizycja określonej ilości obrazów lub akwizycja dokonywana w sposób ciągły (do bufora w pamięci lub na dysku). Ustawienie pobieranej ilości obrazów można dokonać w zakładce General - opcja Frames Per Trigger.
Domyślną wartością jest 1 obraz.

Wyjaśnienia wymaga również pojęcie Trigger. Jest to po prostu informacja o początku akwizycji, inaczej wyzwalanie akwizycji. Sposób wyzwalania można ustawić w zakładce Triggering. Tutaj pozostawiona została wartość domyślna - wyzwalanie natychmiastowe, tzn. równoczesne z początkiem akwizycji (naciśnięciem przycisku Start Acquisition).

Aby móc w pełni korzystać z urządzenia akwizycji, w pewnych sytuacjach konieczne jest ustawienie parametrów tego urządzenia. Przykładowo - wybór źródła sygnału (S-video, composite…). W narzędziu IMAQTOOL można to zrobić korzystając z opcji w zakładce Device Properties.

Na potrzeby akwizycji jednego obrazu pozostałe parametry można pozostawić jako domyślne. Aby sprawdzić czy urządzenie zostało poprawnie skonfigurowane można skorzystać z podglądu (przycisk Start Preview). Jeśli wszystko zostało poprawnie skonfigurowane w oknie Preview powinien pojawić się aktualny obraz z kamery. Uruchomienie procesu akwizycji może chwilę potrwać.

Po zamknięciu podglądu, możliwe jest pobranie obrazu do środowiska MATLAB. W tym celu należy nacisnąć przycisk Start Acquisition i poczekać aż obraz zostanie pobrany (przycisk Start Acquisition stanie się powtórnie aktywny). W przypadku gdy ustawione byłoby wyzwalanie ręczne (Trigger = manual) akwizycję można rozpocząć naciskając przycisk Trigger.
Obraz po akwizycji znajduje się w buforze. Parametry tego bufora można ustawić w zakładce Logging.

Aby przenieść obraz z bufora do przestrzeni roboczej MATLABa, należy nacisnąć przycisk Export Data oraz ustawić nazwę zmiennej, która będzie przechowywała obraz.

W ten sposób osiągnięty został postawiony cel - pobranie obrazu z kamery do przestrzeni roboczej MATLAba. Na zakończenie dodam, że ustawiona w przedstawionym narzędziu konfiguracja może zostać wyeksportowana do późniejszego wykorzystania. Możliwy jest eksport do przestrzeni roboczej MATLABa (jako zmienna) lub do pliku MAT - menu File>Export Hardware Configuration. Możliwa jest również generacja m-kodu realizującego konfigurowanie urządzenia - menu File>Generate M-file.
Posted in MATLAB | 1 komentarz »
czerwiec 9th, 2008 milosz
Zdecydowana większość użytkowników środowiska MATLAB dochodzi w swoje pracy do momentu, w którym musi zaimportować i przetworzyć jakieś dane „z zewnątrz”. Dla większości z nas taka operacja sprowadza się najczęściej do zaimportowania pliku w odpowiednim formacie, z czym MATLAB dzięki zaimplementowanym „wizardom” radzi sobie bardzo dobrze. Nie wszyscy zdają sobie sprawę, że wczytywanie plików to nie jedyna metoda pobierania danych do środowiska. MATLAB potrafi bez problemu podpiąć się do bazy danych lub porozumieć z umieszczoną w komputerze kartą akwizycji danych tudzież podłączonym z zewnątrz oscyloskopem. Pisząc „porozumieć” mam na myśli odczyt danych pomiarowych jak również konfiguracje urządzenia. Taka metoda pobierania danych ma sporo zalet, przy czym największą z nich zdaje się być oszczędność czasu. Chciałem dzisiaj napisać kilka słów na temat Data Acquisition Toolbox, przybornika służącego właśnie do obsługi kart akwizycji danych. Pewną ciekawostką może być fakt, że w celu skorzystania z tego toolboxa wcale nie musimy dysponować kartą pomiarową, a konkretniej, jako kartę akwizycji możemy potraktować zainstalowaną w każdym współczesnym komputerze kartę dźwiękową. Naturalnie karta dźwiękowa, nawet dobra, ma sporo ograniczeń, niemniej jednak:
- Może służyć w celach próbnych oraz do testowania algorytmu, który przy przeniesieniu na komputer z „prawdziwą” kartą akwizycji nie podlega prawie żadnym modyfikacjom.
- Jest wystarczająco dobrym urządzeniem, jeśli chcemy zarejestrować … sygnał dźwiękowy.
Przykład opisany poniżej będzie dotyczył konfiguracji karty dźwiękowej właśnie pod kątem rejestrowania dźwięku.
Pierwszą komendą, jako warto zapamiętać przymierzając się do pracy z Data Acquisition Toolbox jest ‘daqhwinfo’. Wpisana bez żadnych dodatkowych parametrów może dostarczyć nam informacji o zainstalowanych w systemie urządzeniach akwizycji.
-
>> info = daqhwinfo
-
>> info.InstalledAdaptors

Ukazała się lista wszystkich wykrytych w systemie urządzeń akwizycji danych. Nas najbardziej interesuje oczywiście pozycja ‘winsound’, jednak warto zwrócić uwagę również na ‘paraller’. Okazuje się, że jako urządzenie akwizycji możemy traktować również dosyć archaiczny współcześnie port drukarki.
Czas przejść do pracy. Na wstępnie musimy stworzyć w przestrzeni MATLABA obiekt odpowiedzialny za porozumiewanie się z wybranym urządzeniem akwizycji. Wywołanie funkcji daqhwinfo z parametrem ‘winsound’ zwróci nam podpowiedź, jak tego dokonać

Skoro wiemy już jak wygląda konstruktor, piszemy:
-
>> ai = analoginput(‘winsound‘, 0);
Z obiektem musimy teraz skojarzyć kanały, po którym będzie następować akwizycja. Ponieważ mamy możliwość nagrywania na dwa kanały (stereo), skorzystamy z tego.
-
>> chan = addchannel(ai, 1:2);
W zasadzie możemy już zacząć pobierać dane
Operacja sprowadza się do dwóch komend odpowiedzialnych kolejno za rozpoczęcie procesu:
oraz za przeniesienie zarejestrowanych danych do przestrzeni roboczej MATLABA:
-
>> [data, time] = getdata(ai);
Podejrzewam, że niewielu użytkowników zadowoli domyślna konfiguracja karty. Podstawowe informacje o ustawionych parametrach akwizycji możemy uzyskać wpisując w oknie poleceń nazwę utworzonego wcześniej obiektu.
Natomiast bliższe informacje na temat możliwości konfiguracyjnych zdobywamy wykorzystując znane już chyba polecenie

Zaczniemy od zmiany częstości próbkowania zapisywanego sygnału. Wszelkich zmian dokonujemy przy pomocy polecenia ’set’
-
>> SampleRate = 22050;
-
>>
set(ai, ‘
SampleRate‘, SampleRate, ‘
SamplesPerTrigger‘, SampleRate
)
Ustawiając parametr ‘SamplesPerTriger’ pośrednio decydujemy o czasie trwania akwizycji. Ponieważ ilość rejestrowanych próbek jest równa co do wartości częstotliwości próbkowania, proces akwizycji będzie w tym wypadku trwał dokładnie jedną sekundę.
Przejdźmy do nieco bardziej zaawansowanego problemu. Przy bieżących ustawieniach chcąc zarejestrować sygnał musimy w odpowiednim momencie wpisać komendę start(ai). Nie da się ukryć, że jednoczesne pisanie na klawiaturze i śpiewanie do mikrofonu może być to trochę męczące, a w pewnych przypadkach niemożliwe do wykonania. Dlatego też wykorzystamy możliwość zastosowania programowego wyzwalania akwizycji. Ustalimy, że rejestracja ma się rozpocząć automatycznie, w momencie pojawienia się wejściu karty dźwiękowej sygnału (napięcia) o określonej amplitudzie (np. 0,2V).
-
>>
set(ai, ‘
TriggerType‘, ‘
Software‘, ‘
TriggerChannel‘, chan
);
-
>>
set(ai, ‘
TriggerCondition‘, ‘
Rising‘, ‘
TriggerConditionValue‘,
0.2);
W zasadzie musimy uporać się z jeszcze jednym problemem. Przy ustawieniach jak powyżej narastające zbocze sygnału wyzwala rozpoczęcie rejestracji. Nas jednak może interesować również krótki fragment narastającego sygnału, który nie osiągną jeszcze zadanej wartości. Czy jest sposób, by ten fragment również zarejestrować? Tak, jest to możliwe
-
>>
set(ai, ‘
TriggerDelay‘,
-0.5);
Możemy pójść jeszcze krok dalej. Nic nie stoi na przeszkodzie, by wyzwalacz wykorzystany do rozpoczęcia procesu akwizycji był również odpowiedzialny za wywołanie dalszych funkcji, których przeznaczeniem mogłaby być np. analiza rejestrowanych danych.
-
>>
set(ai, ‘
TriggerFcn‘, @daqcallback
);
To tylko przykład, funkcja daqcallback zwraca po prostu czas rozpoczęcia akwizycji. Nic nie stoi jednak na przeszkodzie w zastąpieniu jej funkcją inną, zdefiniowaną przez użytkownika.
To w zasadzie wszystko, co chciałem Wam dzisiaj przekazać. Wspomnę jeszcze tylko, że Data Acquisition Toolbox dodaje również kilka bloczków do Simulinka, umożliwiając pracę z kartami akwizycji (w tym również z kartami dźwiękowymi) także w tym środowisku. Ale to już temat na kolejny odcinek.
Posted in MATLAB | 3 komentarzy »
czerwiec 3rd, 2008 pawel
Ostatnio zostałem poproszony przez użytkownika forum o napisanie bloga na temat układu różniczkowego II rzędu, układ masa – tłumnik – sprężyna. Jest on opisany równaniem różniczkowym F =mx”+ cx’+ kx (apostrof oznacza pochodną), gdzie m - masa, c - stała tłumienia, k - współczynnik sprężystości sprężyny.

Przekształcamy równanie F =mx”+ cx’+ kx do postaci x’’ = 1/m * ( F – cx’ – kx ).
Dlaczego modelujemy układ różniczkujący za pomocą bloków całkujących i przekształcamy to do postaci pokazującej drugą pochodną wyjaśnione zostało w blogu podstawy Simulinka, cz.2 http://www.matlab.matlab.pl/?p=8
Na początku tworzymy nowy model Simulinka i przeciągamy z biblioteki Simulink->Continuous dwa bloki Integrator.

Proszę zwrócić uwagę jak podpisane są sygnały pomiędzy bloczkami, odpowiednio druga pochodna x’’ = przyśpieszenie, pierwsza pochodna x’ = prędkość i wartość x przemieszczenie.
Mamy zamodelowane pochodne, teraz przystępujemy do modelowania stałych 1/m, c i k.Do tego użyjemy bloku Gain z biblioteki Simulink-> Math Operations. Kolejno wypełniamy te bloki stałymi tak jak na modelu poniżej.

Z tej samej biblioteki przeciągamy blok Sum. Dwukrotnie naciskamy na nim i konfigurujemy go tak jak na poniżej.

Następnie podłączamy wyjścia z bloków Gain do Sumatora tak jak poniżej.

Do bloku Sumatora do + (równanie wyżej) podłączamy blok Step z biblioteki Simulink->Sources , do podglądnięcia wartości przemieszczenia masy podłączamy blok Scope z biblioteki Simulink->Sinks i podłączamy tak jak poniżej.

Następnie przechodzimy do Command Window, Matlaba i przypisujemy zmiennym m,k,c odpowiednie wartości.

I uruchamiamy symulacje.

Rozszerzamy czas symulacji do 100.

Zmieniamy solver z ode45 na ode23, aby dokładniej przesymulować ten układ i widzimy różnice na wykresie.

Widzimy jak układ reaguje na wymuszenie skokowe o wartości 1. Parametry m,k,c oraz wartość wymuszenia można zmieniać, otrzymując różne odpowiedzi. Możemy także rozbudować układ o regulator, zmieniać typ wymuszenia itd. Zachęcam do eksperymentowania.
Posted in Simulink | 2 komentarzy »
maj 27th, 2008 pawel
Silnik prądu stałego jest silnikiem elektrycznym zasilanym prądem stałym i służy do zamiany energii elektrycznej na energię mechaniczną.
Jako maszyna elektryczna prądu stałego może pracować zamiennie jako silnik lub prądnica. W tym drugim przypadku wirnik napędzany jest energią mechaniczną dostarczona z zewnątrz, a na zaciskach uzwojenia twornika odbierana jest wytworzona energia elektryczna.
(źródło Wikipedia http://pl.wikipedia.org/wiki/Silnik_pr%C4%85du_sta%C5%82ego)
Do zamodelowania silnika DC użyjemy elementów bibliotek programu Simulink :
1. Część elektryczna w SimElectronics
2. Cześć mechaniczna w SimMechanics.
Biblioteki te razem z SimHydraulics i SimDriveline, tworzą platformę na której możemy modelować systemy fizyczne.

Tworzymy czyste okno Simulinka.

1.Przystępujemy do budowy części elektrycznej. Z biblioteki Simulinka Simscape > Foundation Library > Electrical > Electrical Sources pobieramy blok DC Voltage Source i umieszczamy je na czystym schemacie. Jest to blok generujący sygnał napięciowy o wartości 1 V w obwodzie.
2.Z biblioteki Simscape > Foundation Library > Electrical > Electrical Sensors przeciągamy blok Current Sensor, który przekształca natężenie prądu w sygnał fizyczny proporcjonalny do natężenia.
3.Łączymy ze sobą dwa bloczki, trzymając PPM. Nasz układ wygląda wtedy tak.

4. Mamy źródło prądu stałego, które płynie w obwodzie. W silniku prądu stałego zasilane są magnesy, które indukują pole elektromagnetyczne. Aby uwzględnić rezystancję obwodu oraz indukowane pole, przeciągamy rezystor i cewkę do obwodu, Simscape > Foundation Library > Electrical > Electrical Elements zaznaczamy i przeciągamy na projekt Resistor i Inductor. Układ ten musimy także podłączyć do masy w tym celu przeciągamy także z w/w lokalizacji Electrical Reference.

5. Do układu podłączamy także blok Solver Configuration, PS-Simulink Converter z Simscape > Utilities oraz Scope z Simulink > Commonly Used Blocks Solver Configuration definiuje ustawienia algorytmu do przesymulowania układu elektrycznego. Blok PS-Simulink Converter przekształca sygnał prądowy na sygnał programu Simulink.

6. W tym układzie należy jeszcze wybrać odpowiednio dokładny solver (ustawiamy ode23t stosowany przeważnie do symulowania układów elektrycznych).

7. Po dwukrotnym naciśnięciu na blok PS-Simulink Converter, otworzy się okienko z listą rozwijalną. Z listy wybieramy A , czyli Ampery.

8. Następnie uruchamiamy symulację i na wykresie Scope obserwujemy jej wynik.

9. Zamodelowaliśmy część elektryczną silnika prądu stałego. Teraz dodamy część mechaniczną, która będzie przekształcała energię elektryczną w energię mechaniczną. W tym celu dołączamy do modelu część z biblioteki Simscape > Foundation Library > Electrical > blok Rotational ElectromechanicalConverter. Następnie z biblioteki Simscape > Foundation Library > Mechanical > Rotational Elements bierzemy dwa bloki Inertia, który reprezentuje idealną inercję oraz Mechanical Rotational Reference, blok odniesienia dla układów mechanicznych. Łączymy je ze sobą tak jak poniżej i symulujemy.

Wykres prądowy takiego układu.

Teraz dołączymy jeszcze bloczki, które umożliwią podglądnięcie prędkości wirowania części mechanicznej. W tym celu dołączamy do modelu bloczek z biblioteki Simscape > Foundation Library > Mechanical > Mechanical Sensors > blok Ideal Rotational Motion Sensor. Następnie z modelu kopiujemy bloczki PS-Simulink Converter wraz z wykresem Scope. Dwukrotnie naciskamy na bloczek PS-Simulink Converter i w polu Output Signal Unit wpisujemy rpm. Łączymy ze sobą wszystkie bloki tak jak poniżej jest to przedstawione. Blok referencyjny został skopiowany wyżej dla estetyki.

Następnie symulujemy ten układ i obserwujemy na wykresie do jakiej prędkości rozpędza się część mechaniczna silnika.

Następnie do tego modelu w 4 s dołożymy obciążenie w wysokości -0,1. W tym celu dołączamy do modelu bloczek z biblioteki Simscape > Foundation Library > Mechanical > Mechanical Sensors > blok Ideal Torque Source. Dodajemy także blok z biblioteki Simulink>Sources Step, który w 4 sekundzie wygeneruje nam obciążenie. Konfigurujemy go w ten sposób.

Do modelu elektrycznego chcemy zadać sygnał fizyczny, dlatego musimy go przekształcić. Będziemy potrzebować tak jak poprzednio blok przekształcający, tyle że teraz w drugą stronę. Z biblioteki Simscape > Utilities przeciągamy blok Simulink-PS Converter. Łączymy ze sobą wszystkie bloki tak jak to jest przedstawione niżej.

I symulujemy.

Widzimy jak w 4 s zmienia się sytuacja, prędkość obrotów silnika spada, a natężenie prądu wraca do wartości bliskiej 1.
Model jest gotowy. Teraz możemy dowolnie zmieniać i testować różne parametry silnika, wartości prądu, czy obciążenie, ale to pozostawiam już użytkownikom Simulinka.
Posted in Simulink | 7 komentarzy »
maj 20th, 2008 marcin
Ostatnio pisałem, o tym że Simulink służy do modelowania ukladów dynamicznych. Jedną z metod opisu takiego układu jest wykorzystanie do tego równań różniczkowych. Układy takie nazywamy układami ciągłymi i do ich rozwiązania używamy metod numerycznych.
Jak wygląda równanie różniczkowe każdy widzi. Problem pojawia się gdy musimy przedstawić go za pomocą bloczków. Jest to oczywiście bardzo proste pod warunkiem, że się już wie jak się za to zabrać.
A więc mamy dane równanie w postaci:
-
x’(t) + a*x(t) + b = u(t)
Ten apostrof wyraża kropkę nad x czyli pierwszą pochodną po czasie. Jest to bardzo proste równanie liniowe, które zawsze moglibyśmy przepisać do postaci z macierzą stanu i szybko rozwiązać w ‘czystym’ Matlabie. Różnica polega na tym, że skomplikowane równania nieliniowe w Simulinku rozwiązuje się równie prosto co powyższy przykład.
Teraz najważniejsze: DO MODELOWANIA RÓWNAŃ RÓŻNICZKOWYCH W SIMULINKU UŻYWAMY BLOCZKA CAŁKOWANIA (’Integrator’ w bibliotece ‘Continuous’). Wykorzystywanie do tego bloczka różniczkowania (’Derivative’) jest dużym błędem!!! Pamiętajcie proszę, że Simulink został tak zaprojektowany i tak napisany, aby równania różniczkowe rozwiązywać w nim w taki właśnie sposób jak na rysunku poniżej.

Proszę zwrócić szczególną uwagę na podpisy sygnałów przy bloczku całkującym. Bloczek całkuje więc jeśli na jego wejście wprowadzę różniczkę sygnału x(t) to na wyjściu pojawi się samo x(t). ‘Składanie’ równania rozpoczynam od przepisania go do postaci:
-
x’(t) = - a*x(t) - b + u(t)
Następnie na pusty model przeciągam jeden bloczek ‘Integrator’ i podpisuje jego wejście i wyjście (żeby się nie pomylić). Zgodnie z równaniem pochodna z x po czasie do suma kilku elementów więc podpinam do pochodnej sumator i doprowadzam jego składowe. Efekt jak widać jest bardzo prosty i czytelny, a na dodatek rewelacyjnie się symuluje.
Warunki początkowe równania ustalamy w oknie paramterów bloczka całkowania.

Jak widać na załączonym obrazku, można tam ustawić także wiele innych rzeczy. Bardzo pożyteczne jest ograniczenie na nacałkowanie (”saturation limit’), zewnętrzny reset (’external reset’) czy dodatkowe wejście do bloczka na wartość inicjalizującą (”Initial condition source’ ustawiamy na ‘External’).
A teraz trochę bardziej skomplikowany przykład:
-
x”(t) - u*(1-x(t)^2)*x‘(t) + x(t) = 0
i jego rozwiązanie:

Tym razem mamy do czynienia z ukladem nieliniowym drugiego rzędu, a dokładnie z równaniem oscylatora Van der Pola.
Zasada działania ta sama:
- ustalić rząd równania i przeciągnąć tyle bloczków całkujących co rząd
- podpisać sygnały wchodzące i wychodzące z bloczków całkujących
- jeśli jest taka możliwość to uzyskać pochodne niższego rzędu z pochodnych wyższego rzędu np. prędkość z przyspieszenia itp.
- przepisać równanie tak aby po lewej stronie pozstała tylko pochodna nawyższego rzędu
- przy pomocy sumatora i innych niezbędnych operacji konstruujemy prawą stronę równania
Na koniec należy jeszcze zdefiniować odpowiedni Solver i możemy oglądać piękne charakterystyki. Tutaj drobna uwaga: o doborze solvera będziemy jeszcze pisać, ale póki co pamiętajcie, że układy ciągłe (opisane równaniami różniczkowymi) symuluje się przy pomocy sloverów zmienno-krokowych.
Posted in Simulink | 3 komentarzy »
maj 6th, 2008 marcin
Żeby rozwiązać wszelkie wątpliowści i nieporozumienia, które pojawiają się wśród użytkowników naszego ulubionego środowiska postanowiłem napisać do czego tak naprawdę służy Simulink.
Definicja mówi, że Simulink jest to narzędzie służące do symulacji układów dynamicznych. Te “układy dynamiczne” układa się z bloczków i w nomenklaturze Simulinka nazywa “modelem”. Modele możemy symulować tzn. wyliczać wartości stanów i wyjść w czasie. Za proces symulacji odpowiada narzędzie nazywane “solverem”. To on odpowiada za metody numeryczne wykorzystywane w trakcie symulacji i jej parametry (krok, precyzja, dopuszczalne błędy itp.).
Przejdźmy teraz do bardziej szczegółowych definicji.
Model - zestaw bloczków opisujących układ dynamiczny, gdzie przez układ dynamiczny rozumiemy coś co można opisać przy pomocy trzech typów równań: równania algebraicznego, równania stanu w wersji ciągłej (równanie różniczkowe) i równania stanu w wersji dyskretnej (równanie różnicowe).
Te trzy typy równań można ze sobą oczywiście dowolnie mieszać tzn. na jednym modelu możemy umieścić równocześnie równanie różniczkowe i różnicowe opisujące inne fragmenty ukadu.
Projektanci The Mathworks projektując Simulinka zakładali, że tworzenie modelu może przebiegać na dwa sposoby:
- modelowanie algorytmu z wyraźnie widocznym strumieniem (kierunkiem), w którym przesuwają się dane od wejścia do wyjścia, typowe dla algortymów przetwarzania sygnałów, telekomunikacji itp.
- modelowania układu sterowania tzn. modelujemy obiekt rzeczywisty (najczęściej ciągły, a więc używamy równań różniczkowych) i regulator (najczęściej dyskretny, a więc używamy równań różnicowych).
O tym jak modelować poszczególne typy równań zostanie omówione w kolejnych wpisach na blogu, ponieważ nie zawsze jest to trywialne. Teraz wymienie jeszcze kilka cech modeli, które być może skuszą Was do zainteresowania sie tym narzędziem:
- Modelowane równania różniczkowe mogą być równaniami nieliniowymi i jest to jedna z ważniejszych cech tego programu. Właściwie zaimplementowanie równania nieliniowego, a liniowego w Simulinku niczym się nie różni, metodyka postępowania jest taka sama. Myślę, że to jest podstawowy powód, który może przekonać fantyków MATLABa, aby zajrzeli również do Simulinka.
- Podczas symulacji równań dyskretnych możemy ze sobą mieszać elementy o różnych czasach dyskretyzacji, co w efekcie pozwala nam zasymulować zaawansowany wieloprocesowy system sterowania przed jego zaimplementowaniem np. na mikrokontrolerze lub w systemie czasu rzeczywistego.
- Jeśli nieobca Ci jest notacja maszyn skończenie stanowych na pewno się ucieszysz, że Stateflow pozwali Ci ją wykorzystać.
- Modele SImulinka mogą wysyłać dane do workspace MATLABa jak i je stamtąd pobierać. W ogóle sama symulacja może być wywoływana z poziomu skryptu, tak jak jej parametry mogą być ustawione przez skrypt. Całość operacji jest wykonywana w tle, bez otwierania Simulinka więc zaawansowany programista MATLABa może potraktować modele jako funkcje o dziwacznym interfejsie.
- I ostatnie - tak to prawda - z modeli Simulinka można wygenerować kod w języku C, a nawet HDL (VHDL/Verilog) używając do tego odpowiednich rozszerzeń tzn. Real-Time Workshop i Simulink HDL Coder. I od razu rozwieję wątpliwości: kod ten jest bardzo czytelny i bardzo łatwo integruje się go z własnymi projektami.
Mam nadzieję, że zainteresowałem wszystkich Simulinkiem. W następnej części omówię jak składać równania różniczkowe.
Posted in Simulink | 2 komentarzy »
kwiecień 1st, 2008 Marcin Piątek
Często początkujący użytkownicy Matlaba mają problem ze zrozumieniem jak działają niestandardowe opartatory pojawiające się w Matlabie. Chodzi o operator zakresu “:” (jest to moja autorska nazwa więc proszę się nie oburzać) i opartor operacji skalarnych “.” (podobna historia).
Zacznijmy jednak od początku… Podstawowym typem zmiennych w Matlabie jest: macierz (matrix). Gdy tworzymy nową zmienną w edytorze lub linii poleceń to możemy być pewni, że zostanie stworzona właśnie macierz. Nawet jeśli podamy tylko jedną skalarną liczbę jako wartość inicjalizującą, to tak naprawdę powstanie macierz o rozmiarach 1×1.
-
>> A = 2
-
-
A =
-
-
2
-
-
>> A = [1 2; 3 4]
-
-
A =
-
-
1 2
-
3 4
Zastosowana precyzja obliczeniowa każdego z elementów tej macierzy to double. Użytkownik może jednak wymusić zastosowanie innego typu poprzez funkcję rzutującą.
Jeśli teraz chcemy stowrzyć np. tablicę składającą się z 10 kolejnych elementów od 1 do 10 to możemy to wykonać na dwa sposoby: po pierwsze wymieniejąc te elementy po drugie korzystając z operatora zakresu.
-
>> A = 1:10
-
-
A =
-
-
1 2 3 4 5 6 7 8 9 10
-
-
>> A = 1:2:20
-
-
A =
-
-
1 3 5 7 9 11 13 15 17 19
Jak widać na powyższym przykładzie można stosować także ten operator podając krok (domyślnie zakładamy 1).
Ten sam operator możemy wykorzystać do wyciągania elementów z macierzy np.
-
-
-
A =
-
-
17 24 1 8 15
-
23 5 7 14 16
-
4 6 13 20 22
-
10 12 19 21 3
-
11 18 25 2 9
-
-
>> A(1,2)
-
-
-
-
24
-
-
>> A(1:3,2)
-
-
-
-
24
-
5
-
6
Należy pamiętać o tym, że wyciągając poszczególne elementy macierzy przy pomocy nawiasów okrągłych piszemy NAJPIERW WIERSZ, A POTEM KOLUMNĘ. No i numerujemy od 1 nie od 0. Ostatnie polecenie w porzednim przykładzie mozna oczywiście zastąpić wersją:
ale nie wygląda to ładnie, a już na pewno nie jest praktyczne.
Pozostaje jeszcze sytuacja gdy chcemy wyciągnąć cały wiersz lub kolumnę. Spróbujmy np. wyciągnąć cały drugi wiersz poprzednio zdefiniowanej macierzy:
-
>> A(2,1:5)
-
-
-
-
23 5 7 14 16
-
-
>> A(2,:)
-
-
-
-
23 5 7 14 16
Ta druga wersja jest bardzo wygodna i często stosowana, ale sprawia początkującym najwięcej problemów więc przypomnę: najpierw podajemy wiersz (w naszym przypadku 2), a potem kolumnę (w naszym przypadku elementy z zakresu: pierwszy element do ostatani element). Analogicznie wyciągniemy całą kolumnę:
lub całą macierz (przykład BARDZO niepraktyczny)
-
>> A(:,:)
-
-
-
-
17 24 1 8 15
-
23 5 7 14 16
-
4 6 13 20 22
-
10 12 19 21 3
-
11 18 25 2 9
Pozostaje jeszcze kwestia kropki “.”. Oparator ten służy do wykonywania operacji “tablicowych, skalarnych?” na macierzach i tablicach. Wyobraźmy sobie sytuację, że mamy podane dwa wektory:
-
>> A = [1:10]
-
-
A =
-
-
1 2 3 4 5 6 7 8 9 10
-
-
>> B = [1:2:20]
-
-
B =
-
-
1 3 5 7 9 11 13 15 17 19
i potrzebujemy pomnożyć każdy element wektora A razy każdy element wektora B. I do tego służy właśnie ten operator tzn.
-
>> A.*B
-
-
-
-
1 6 15 28 45 66 91 120 153 190
To samo tyczy sie np. dzielenia, ale najbardziej pożyteczne jest przy potęgowaniu, bo gdy napiszemy:
-
-
-
A =
-
-
17 24 1 8 15
-
23 5 7 14 16
-
4 6 13 20 22
-
10 12 19 21 3
-
11 18 25 2 9
-
-
>> A^2
-
-
-
-
1090 900 725 690 820
-
850 1075 815 720 765
-
700 840 1145 840 700
-
765 720 815 1075 850
-
820 690 725 900 1090
to policzymy kwadrat macierzy. Jest to operacja co prawda bardzo przydatna i fajnie, że Matlab ją potrafi wykonać, ale znacznie częściej będziemy potrzebowali podnieść po prostu do kwadratu poszczególne elemnty macierzy (jeśli ktoś nie rozumie na czym polega różnica odsyłam do podręcznika z algebry liniowej).
-
>> A.^2
-
-
-
-
289 576 1 64 225
-
529 25 49 196 256
-
16 36 169 400 484
-
100 144 361 441 9
-
121 324 625 4 81
Posted in MATLAB | Brak komentarzy »