1 Mandelbrot

a) Skriv en funktion mandelbrot(c) som tar det komplexa talet c och beräknar följande:

låt z börja med värdet 0

returnera true om nedanstående beräkning kan upprepas 100 gånger utan att absolutbeloppet av z blir större än 2

z = z*z + c

b) Skapa en matris M med nollor, av storleken 401 * 401, med hjälp av numpy.zeros()

låt a ta alla värden mellan -2 och 2, med stegstorleken 0.01.

låt b ta alla värden mellan -2 och 2, med stegstorleken 0.01.

för alla kombinationer av värden på a och b, beräkna mandelbrot(a + bi). Om värdet blir true så lagra värdet 1 i matrisen M på den position som har indexet [(a+2)*100, (b+2)*100], annars får nollan stå kvar.

Titta på hur M ser ut. Man kan visa matrisen som en svartvit bild med hjälp av matplotlib.pyplot.imshow(M, cmap='gray', extent=(amin, amax, bmin, bmax)). 'extent' används för att ge rätt skala på axlarna - amin ska vara det minsta värdet på a, amax det största, osv, Varje pixel motsvarar nu värdet i den positionen, antingen svart eller vit.

(Frivilligt tips: i stället för att lagra ettor eller nollor så kan ni modifiera madelbrot(c), och lagra hur många steg ni tog utan att beloppet av z överstiger 2, då kan ni göra det till en lite vackrare fraktal, genom att ansätta en fin färgskala. Tips finns i tutorial till Pyplot).

c) Välj ut ett nytt del-intervall av a-värden och b-värden, sådant att

skillnaden mellan största och minsta värde på a resp b är en tiondel så stort som förra gången.
Intervallen fortfarande är indelade i 400 steg.
Det nya intervallet innehåller både vita och svarta pixlar från föregående steg

(Ni kan se det som att ni plockar ut en liten ruta ur den föregående bilden, och förstorar up den.)

Upprepa beräkningen på det nya intervallet, och plotta hur det ser ut.

Upprepa ett par gånger till, med nya intervall, som spänner en tiondel av föregående intervall.

OBS:
Om ni väljer de nya intervallen enligt den röda rutan nedan, så blir den inzoomade bilden tråkig, eftersom den bara innehåller pixlar av samma färg:

Om ni i stället väljer ett intervall som innehåller både svarta och vita pixlar, så kommer den inzoomade bilden att vara mer intressant att titta på:

Vilka värden på a resp. b som ger en intressant ruta att zooma in i, kan ni se i den första bilden ni plottar.

För att få rätt skala på bilden (så att det inte bara står pixelindex 1-401 längs med axlarna), kan ni använda er av argumentet extent i imshow, såhär:

plt.imshow(M, cmap='gray', extent=(amin, amax, bmin, bmax))

Där amin är minsta värdet på a, amax största värdet på a, o.s.v.

2 Bildeffekter

I den här uppgiften ska ni använda 2-dimensionell faltning, och omvandla en bild till en "tuschteckning".

Detta består av följande steg

Konvertera bilden från färg till svartvitt.
Använd faltning för att hitta var derivatan i bilden är som störst, dvs, var det är störst skillnad mellan intill-liggande pixlar. Dessa områden ska bli strecken i "teckningen".
Lägg till nya pixlar runt de som är kvar, för att få tjockare streck.

Nedan följer stegen i mer detalj:

1) Ladda in en bild med hjälp av scipy.ndimage.imread (i senaste versionen av SciPy har denna funktion flyttats till matplotlib.pyplot.imread). Ni kan använda vilken bild som helst, t.ex ett foto på er själva. I dessa instruktioner använder vi den här bilden på KTH, men ni måste visa att er kod fungerar på en annan bild än denna. Skapa en ny ndarray som innehåller medelvärdet av de tre grundfärgerna i ursprungsbilden, så att det blir en svartvit bild: Borggårdedn i svart-vitt

Tips: om ni har svårt att få en snygg svart-vit bild med pyplot.imshow(), så kan ni behöva ange specifikt att ni vill visa bilden som en gråskala. Hur man gör det hittar ni i dokumentationen till imshow.

2) Falta bilden med de sk. Sobel-kernels, som beräknar derivatan i x- respektive y-led.

Gx = , Gy =

Spara två separata ndarray:er, en som är bilden faltad med Gx, en som är bilden faltad med Gy, enligt följande modell, där derivatorna i x- rep. y-led sparas i xd resp. yd:

xd = ndimage.convolve(bild_svartvit, Gx, mode='constant')
yd = ndimage.convolve(bild_svartvit, Gy, mode='constant')

För att beräkna var våra streck ska vara så bryr vi oss inte om i vilken riktning som derivatan går, dvs om det blir mörkare uppåt eller neråt eller höger eller vänster, bara hur stor skillnaden är mellan pixlarna just här. Alltså beräknar vi beloppet av derivatan, så att vår linjebild S i varje pixel (i,j) blir $s_{i, j} = \sqrt{x d_{i, j}^{2} + y d_{i, j}^{2}}$ .
Om ni nu plottar S så borde den se ut så här:

3) För att få svarta linjer mot en vit bakgrund kan man gå igenom alla pixlar i S och låta alla pixlar med ett värde över t.ex 100 ersättas med 0, och de under 100 med 255. Då får man en bild som ser ut såhär (beroende på hur bilden ser ut kan andra värden än 100 ge bättre resultat):

Tips: Ni kan ändra på alla element i en array S som uppfyller ett visst villkor genom att skriva såhär:

index = <villkor>
S[index] = <nytt värde>

<villkoret> kan här t.ex vara S>100, och <nytt värde> 0

Det blir mer ännu mer kompakt att skriva:

S[<villkor>] = <nytt värde>, alltså t.ex:

S[S>100] = 0

Som nästa steg kan man låta bilden ritas med en tjockare "penna". Om man t.ex faltar med en nxn-matris där alla element är 1, och sätter alla element där resultatet är 255*n^2 till 255, och övriga till noll, så får man följande bild (n=5):

Prova gärna med olika stora matriser, eller olika värden på de olika konstanterna som ingår!

3. Ljudfiler och fouriertransformer:

I den här uppgiften ska vi titta på ljudinspelningar, och använda några olika verktyg för att analysera och manipulera dem.

a) Läs in och visualisera en ljudfil

Använd er av scipy.io.wavfile.read för att läsa in en fil med en sampling av en pianoton. Plotta ljudvågen. Använd scipy.fftpack.fft för att beräkna fouriertransformen F av ljudsamplingen. Plotta absolutbeloppet av första halvan av F. Hur ska x-axeln skalas för att frekvenserna ska stämma? (tips: titta på samplingens längd och utgå från vilken frekvens som behövs för att så många vågor ska behövas för att bli den tiden. Det i:te elemetet i F motsvarar den frekvens där det ryms i hela vågor under samplingens tid. Om t.ex i=10, och samplingen är 2 sekunder lång, så ska tio vågor bli 2 sekunder, och varje våg måste vara 0.2 s, och frekvensen är alltså 5 Hz).

Sampling att testa. Detta är ett C på ett piano, ca 3 sekunder långt:
Piano_1_C_.wav

b) Läs in en ljudfil och hitta tonen

Skriv en funktion som kan ta en frekvens som indata, och returnera namnet på den ton som ligger närmast, som en sträng. Som exempel ska argumentet 440 returnera "A". Eftersom ni inte behöver kunna ange rätt oktav, är ett tips att dubblera eller halvera värdet tills det ligger inom ett givet intervall, t.ex 440-880, och sedan beräkna vilken ton det är.

Skriv nu ett program som kan läsa in en ljudfil, hitta vilken frekvens som har högst amplitud (absolutbeloppet i fouriertransformen), och använder funktionen ovan för att ange vad tonen heter. Testa på ljudfilerna nedan. För att kontrollera att ert program fungerar kan ni anta att den första filen innehåller tonen C.

Piano_1
Piano_2
Piano_3
Piano_4
Piano_5
Piano_6

c) Omvandla ett ackord från dur till moll

Nu ska vi ta en mer avancerad sampling, med ett C-durackord (består av tonerna C-E-G), och byta frekvens på enbart mitt-tonen, från E till Ess, så att det blir ett C-mollackord (består alltså av tonerna C-Ess-G). Vi vill alltså inte byta frekvens på C:et eller G:et, så det räcker inte att bara byta till en lägre samplingsfrekvens, och därmed minska tonhöjden på hela ackordet. Vi har förenklat problemet lite, så ha inte för höga förhoppningar på ljudkvaliteten.

Skriv ett program som:

Läser in samplingen cdur.wav.
Visualiserar de ingående frekvenserna, med hjälp av fourier-transformen. Hitta de toppar i frekvensplotten som motsvarar E:et, och dess övertoner (dvs hela multiplar av basfrekvensen). Vi struntar i de frekvenser som ger pianot dess klangbild, och tittar bara på hela oktavsteg.
Flytta värdena i den fouriertransformerade vektorn från E till Ess (Eb). Se till att inte bara få med den absoluta toppen, utan även de kringliggande värden som hör till samma topp. viktigt: Kom ihåg att fourier-vektorn är speglad runt mittvärdet, så ni måste flytta på motsvarande sätt i andra halvan. Om ni har en fouriertransform som är n element lång, och ni hittar toppen på element m, m<n/2, och tänker flytta från m till m-4, så måste ni samtidigt flytta värdet i element n-m till n-m+4:
Plotta frekvenserna för både dur- och mollackordet i samma plot, och kontrollera att det bara är rätt frekvenstoppar som har flyttats. Det borde då kunna se ungefär så här (Den blå plotten är mollackordet, orange plot är dur. Den blå plotten syns bara där den är olik den orangea):
Använd scipy.fftpack.ifft för att invertera fouriertransformen och skapa en ljudsignal utifrån era modifierade frekvenser. På grund av avrundningsfel finns en risk att er ljudsignal innehåller imaginära komponenter, så se till att plocka bort dessa och bara behålla realdelen av varje element. Skriv sedan ljudsignalen till en wav-fil med hjälp av scipy.io.wavfile.write.
Notera att er ljudsignal nu består av floats. Funktionen wavfile.write förväntar sig att floats ska vara på intervallet [-1, 1], så ni kan behöva skala ner dem för att inte få distortion. Ni kan också styra datatypen på wav-filen genom att ange datatyp på array:en med ljudsignalen, med hjälp av ndarray.astype(), t.ex så här:
"data.astype(np.float32)" eller "data.astype(np.int16)".

Nu borde ni kunna spela upp ert c-mollackord. Jämför med durackordet.

Ni kan testa på andra ljudfiler också, om ni vill. Kanske kan ni spela in egna?

Felsökningstips:
Om din wav-fil inte låter när du spelar upp den, kan du testa med ett lite enklare program. Om det enklare programmet inte heller fungerar så kan du testa genom att byta datatyp, och skala om värdena. När det fungerar för sinusvågen så kan du se till att använda samma inställningar för ditt ackord.

Tänk på att wav-filer med float förväntar sig värden mellan -1 och 1, medan int16-versionen vill ha heltal mellan-32768 och 32767. Du hittar värden för alla datatyper här. Om du använder värden mellan -1 och 1 för en int-fil, så kommer din volym bli så låg så att den inte hörs, medan för stora värden i en float-fil ger en massa brus.