北大心理学系Psychtoolbox1

1、PsychtooboxLihan Chen2015.4.28MATLAB 图像刺激 图像处理 图像基础 图像的变换与操纵 计算生成图像图像基础图像基础像素、颜色、索引图像图像基础像素 数字图像 离散的坐标：像素 离散的强度：0255 (黑白)图像基础颜色 RGB颜色系统 每个像素的红、绿、蓝均由uint8表示 (0255) alpha通道：透明度 24位/32位真彩色imtool( photo.jpg);用MATLAB存取并显示图片读取并显示图片 im = imread(tmp.png); size(im) imshow(im)写入图片 imwrite(im, tmp.jpg);a=ones(

2、400,40),zeros(400,40); %先生成一根白条和一根黑条a1=repmat(a,1,5); % 将这个组合复制5遍a2=xor(a1,a1);%将这个黑白条纹与自己的转置做异或逻辑运算imshow(a2) %棋盘格就生成了实例：红蓝3D与双眼竞争 im1 = imread(Koala.jpg); im2 = imread(Penguins.jpg);% 用第二幅图替代第一幅图的红色通道 im1(:,:,1) = im2(:,:,1); % 把绿色通道清零 im1(:,:,2) = 0; imshow(im1)image与colormap image函数是MATLAB最基本的图像

3、显示函数，可以绘制索引图像，即每个像素的值对应颜色查找表中的索引 colormap：定义图像显示用的颜色查找表 imagesc将数据scale后绘制成图（例如绘制相关矩阵）image与imshow Imshow仅用于显示由RGB或灰度值定义的图像（image也可以） 无论是哪个函数，若图像是以uint8表示的，则取值范围为0255，若以double表示，则取值范围是01image 绘制数据 246810123456789100.50.550.60.650.70.750.80.850.90.951imagesc & colormap 绘制相关矩阵a = rand(10);for i = 2:10

4、a(:,i) = a(:,i-1)+rand(10,1)*0.5;endimagesc(corr(a)colorbarcolormap autumn图像的变换与处理图像的变换与处理调整曲线、窗口化像素亮度的变换操作 (PS中的“曲线”)实例：去除扫描时背面透过来的影 im = imread(book-000123.png); g = rgb2gray(im); d = double(g); y = uint8(d*1.7-85); imshow(y); imwrite(y, out.png);图像的“窗口”操纵每个像素的亮度：乘以一个系数生成高斯窗口窗口化图像 imagesc(GaWindow

5、) tmp = double(im).*GaWindow; Windowed = uint8(tmp); imshow(Windowed)邻居处理邻居处理高斯模糊、边缘检测空间滤波 定义一个矩阵(称为滤波器filter) 将该矩阵覆盖在每个像素点上 对所有被覆盖的点 将原图像矩阵和该矩阵的对应点相乘 求和变成新图像的点 由邻居决定新的像素点的强度高斯模糊 改变窗口大小，令x = y = 10, stdev = 5，生成一个“高斯窗口” 让窗口矩阵的和等于1GaFilter = GaWindow/sum(sum(GaWindow); 用imfilter或filter2或conv2对原图进行滤波i

6、mf = imfilter(im, GaFilter);imshow(imf)高斯模糊的图片边缘检测原始图像原始图像 grayedge(gray)imfilter(gray, fspecial(sobel)更多的滤波器 用fspecial生成二维的特殊滤波器 average 简单平均 gaussian 高斯滤波 motion 模拟摄像机抖动 sobel 检测边缘算子 doc fspecial计算生成图像计算生成图像正弦光栅, Gabor Patch, 生成正弦光栅 X, Y=meshgrid(1:50,1:100); imshow(sin(X) imshow(sin(X+Y) imshow(s

7、in(X+Y)/2) imshow(sin(X+2*Y)/2) imshow(0.5+0.5*sin(3*X+4*Y)/6)完整代码x,y=meshgrid(1:Boundx, 1:Boundy); % define the costant for the cosine wo=SAngle/180*pi; % conversion from deg to radPhRad=PhGrad/180*pi;% conversion from deg to radf=Fcxi/xysize(1); % frequency conversion ax = f*cos(wo); % convert the

8、 spatial frequency along yby = f*sin(wo); % convert the spatial frequency along x % Calculate the Sinusoid %(+1 is to obtatin all values being more than 1)ImOut=sin(2*pi*(ax*x+by*y)+PhRad)+1;生成Gabor% 先生成正弦光栅 X,Y=meshgrid(-50:50,-50:50); im = sin(X+2*Y)/3); imshow(im)% 生成窗口“盖住”光栅 Ga = exp(-(X/20).2)

9、- (Y/20).2); imshow(Ga.*im) Psychtoolbox函数介绍（图像刺激的实现） SCREEN函数与视觉刺激函数与视觉刺激产生视觉刺激的三种方式 导入图片(jpg,. tif,) 在PTB里画图 文本图片（text figures)Screen函数 Screen是PTB与视频显示有关的函数的集合 都有哪些子函数？ Command window输入Screen Screen(函数字符串，参数列表） /Screen windowPtr, rect=Screen(Openwindow,windowPtrOrScreenN

10、umber, , color, rect, pixelSize,numberOfBuffers ,stereomode , multisample , imagingmode , specialFlags , , clientRect);使用Screen函数 常用的子函数、使用顺序及功能 查询有多少屏幕 Screens 开窗口（创建主页面） OpenWindow 绘制视觉刺激 画线 DrawLine, 画点 DrawDots, 贴图 DrawTexture 页面切换 Flip 关闭 CloseAll 附注 ShowCursor; HideCursor;Screen函数：Screens子函数 S

11、creens子函数：查询有多少屏幕 FAQ: 有多少屏幕？ 若有多个显示器，一般主屏幕用来控制，编号最大的屏幕用来呈现刺激 此步可省略，可以直接开0号窗口 screens=Screen(Screens); screenNumber=max(screens);Screen函数： OpenWindow子函数 OpenWindow用来打开一个窗口(“画布”) w, rect = Screen(OpenWindow, screenNumber); 返回值 w: 该窗口的“句柄值”，类似于“指针” 以后对这个窗口进行任何操作都要用这个指针 rect: 该窗口的“矩形”，由左上、右下两个点的座标定义，比如

12、0 0 1024 768 更多参数 背景颜色、窗口大小？ Screen OpenWindow?绘制视觉刺激 画什么？ 线条？DrawLine 点？DrawDots 多边形？FillPoly / FramePoly 椭圆？FillOval / FrameOval 文字？DrawText 照片？MakeTexture + DrawTexture 效果？ 位置？ 透明？Screen函数： DrawDots子函数 画点：DrawDots Screen(DrawDots, w, xy, size, color, center,type); 用来画一个或一堆点，它的参数们： w: 窗口的“指针” xy:

13、两行的行向量，每个点的x、y座标各占一行 size: 点的宽度，以像素计 color: (每个)点的颜色( R G B 或 R G B A ) center: 点座标的相对中心 type: 方点，圆点，还是抗锯齿的圆点Screen函数： DrawLine子函数 画线：DrawLine Screen(DrawLine, w, color, fromH, fromV, toH, toV, penWidth); 用来画一条线，它的参数们： w: 窗口的“指针” color: 线的颜色 from/to H/V: 起点/终点的x/y座标 penWidth: 线条的粗细Screen函数： MakeText

14、ure子函数 Texture -由图片转换而来的“纹理”，便于快速贴图 怎么把图片转成纹理？ image = imread(photo.jpg); textureIndex=Screen(MakeTexture, w, image); 返回值textureIndex就是这个纹理的“指针” DrawTexture子函数 把纹理（即照片）画到屏幕上 Screen(DrawTexture, w, textureIndex); 把照片画到哪里？ 源矩阵 src_rect=Screen (Rect,tex); 画纹理的哪部分（默认全画）？ 摆在中心 ctr_rect=CenterRect(src_rec

15、t,wRect); 偏移 dst_rect=OffsetRect(dst_rect,-150,150); 画到偏移完的目标上 Screen(DrawTexture, w, tex, src_rect, dst_rect);Screen函数： Flip子函数 Psychtoolbox有两个缓冲区(“画板”) 屏幕上显示的只是其中一面 显示的同时可以在下面的画板中继续画 需要更新的时候，直接切换，速度快 把画板切换上去 马上就翻：Screen(Flip,w); T时刻(或之后)再翻：vbl = Screen(Flip,w, T); 返回值vbl是真正刷屏时刻的时间戳(GetSecs)Screen函

16、数： CloseAll子函数 为了避免程序运行时出错，而无法执行到CloseAll就卡在中间关不了全屏窗口 try screens=Screen(Screens); screenNumber=max(screens); w, rect = Screen(OpenWindow, screenNumber); Screen (CloseAll); catch Screen (CloseAll); rethrow(lasterror); end呈现透明刺激 Alpha通道 R G B A中第四个通道 0为完全不透明，255为完全透明 设置Blend Function Screen(BlendFunc

17、tion, w, GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); “去掉图片中不要的部分”？Psychtoolbox-声音刺激与反应录入 声音基础与产生声音声音基础与产生声音采样率、位数声音基础 音频信号的表示 1或2行(列)采样得到的行(列)向量 时间上的离散：采样率 存储上的离散：位数 (在MATLAB中是double) 声道：单声道和立体声 采样率：每秒离散信号的采样个数 44100Hz 位数：用多少位的数字表示信号对正弦波的4位采样用MATLAB存取并播放声音读取并播放声音 wav, sf, nbits = wavread(sound.wav); size

18、(wav) t = size(wav,1)/sf % 时长 = 采样点个数/采样率 wavplay(wav, sf)写入声音 wavwrite(wav, sf, nbits, sound.wav);播放声音：wavplay 改变采样率 改变播放时长，并变调 wavplay(wav,8000) wavplay(wav,80000) “异步”播放 指定async参数可以使声音在“后台”播放 否则wavplay会等待播完，然后才执行下一句 利用async可以“叠加”两个声音 wavplay(wav, 44100, async)wavplay 与 sound 类似imshow和image，MATLAB

19、中这两个函数都可以播放一段声音，用法类似 只有wavplay可以“异步”播放 soundsc将数据scale到-1 +1区间后播放（如果有意更改声音的音量，不要用这个函数）生成噪音 生成一段单声道噪音就是生成一列随机数 sf = 44100; % 采样率 t = 1.5; % 持续时长 noise = rand(sf*t, 1); wavplay(noise, sf)plot(noise(1:200)02040608010012014016018020001生成纯音 生成一段单声道纯音就是生成一列正弦波 sf = 44100; % 采样率

20、t = 1.5; % 持续时长 f = 261.6; % 纯音频率(C4 = 261.6Hz) tmp = linspace(0, t, sf*t); % 计算正弦用的采样矩阵 tone = sin(2 * pi * f * tmp); wavplay(tone, sf)plot(tone(1:200)020406080100120140160180200-1-0.8-0.6-0.4-0.60.81生成复合音叠加若干纯音sf = 44100; % 采样率t = 1.5; % 持续时长f = 261.6; % 基准频率(C4 = 261.6Hz)tmp = linspace(

21、0, t, sf*t);tone1 = sin(2 * pi * f * tmp);tone2 = sin(2 * pi * 2*f * tmp);tone3 = sin(2 * pi * 4*f * tmp);harmonic = tone1 + tone2 + tone3; % 叠加并归一化harmonic = harmonic / max(abs(harmonic);wavplay(harmonic, sf)plot(harmonic(1:200)020406080100120140160180200-1-0.8-0.6-0.4-0.60.81生成一段声音拼接若干纯音

22、 sf = 44100; % 采样率 t = 0.5; % 持续时长 f1 = 261.6; f2 = 293.6; f3 = 329.6; % (C4, D4, E4) tmp = linspace(0, t, sf*t); do = sin(2 * pi * f1 * tmp); re = sin(2 * pi * f2 * tmp); mi = sin(2 * pi * f3 * tmp); silence = zeros(1, sf*t); % 静音 wavplay(do re mi silence do re mi re do, sf)声音处理的几个一般函数 纯音: MakeBeep

23、beep, samplingRate1= MakeBeep(freq,duration, samplingRate) freq-音调 duration-音频长度（秒） samplingRate- 采样频率Beep-音频数据返回值 samplingRate1- 返回采样频率 Beeper 格式 Beeper (frequency, loudness, duration) loudness “01” Snd(command, signal,rate,samplesize) command集成- Open, Play, Close, Wait, Quiet, DefaultRate, IsPlayi

24、ng signal-单、双通道 rate-播放频率 samplesize-比特率(8或16） Snd(Play, sin(0:5000) Snd(Close) freq = 44100;tone(1,:) = 0.9 * MakeBeep(1000, 0.04, 44100); tone(2,:) = tone(1,:);toneb(1,:) = 0.9 * MakeBeep(500, 0.04, 44100); toneb(2,:) = toneb(1,:); Snd(Play, tone, freq, 16); Snd(Quiet)Snd(Close) 音频数据播放: sound Soun

25、d(y, Fs) Sound(y,Fs,bits) Fs-采样频率 bits- 比特数 音频播放器：audioplayer player= audioplayer(Y,Fs,nBits) player= audioplayer(recorder, ID) 创建录音机对象：audiorecorder recorder=audiorecorder(Fs, nBits, nChannels, ID)录制一段声音录制一段声音recorder=audiorecorder(44100, 8, 1, 0);disp(Start Recording);Recordblocking(recorder, 6);

26、% record 6 secsdisp(End of Recording);play(recorder);audiodata=getaudiodata(recorder);plot(audiodata);产生其他波形 chirp: 产生信号的频率随时间的增长而变化 t=0:0.005:2; % 2 secs 5000 Hz SF y=chirp(t,220,2,440); wavplay(y,5000) sawtooth: 锯齿波 fliplr, flipud: 把矩阵左右或上下颠倒 repmat: 重复矩阵频谱图“通过按键音还原手机号码” 频谱图 t=0:0.0001:2; % 2 secs

27、 10000Hz y=chirp(t,220,2,4400); wavplay(y,10000) spectrogram(y,128,120,128,1e5)立体声、ITD和ILD 立体声由两列-1/+1之间的数表示 ITD(两耳时间差)和ILD(两耳强度差)是人进行声源定位的两个线索 HRTF(头相关传输函数)刻画人耳接收空间中某一点传来的声音的特征 /projects/HRIRrepository声音的操纵声音的操纵包络与滤波音量控制 数字音频中，0dB是最大的音量，在此基础上可以减小音量，如-10dB wav * 10(db/20) 可将

28、wav音量降低为db分贝，db0声音的淡入和淡出 在声音前后一小段加上“门控” gatedur = .01; gate = cos(linspace(pi, 2*pi, fs*gatedur); gate = (gate + 1) / 2; offsetgate = fliplr(gate); sustain = ones(1, (length(tone)-2*length(gate); envelope = gate, sustain, offsetgate; smoothed_tone = envelope .* tone;声音的包络(envelope) 获取音频设备信息: audiode

29、vinfo devinfo= audiodevinfo Name= audiodevinfo(IO, ID); ID= audiodevinfo(IO,rate, bits, chans)Low latency, precisely timed sound with PsychPortAudio Based on PortAudio, a free, open-source, cross-platform audio library for realtime audio: http:/ Features: Multi-channel playback, recording and full-d

30、uplex feedback operation Low-latency sound output on MacOS/X, and MS-Windows (ASIO- required) Precisely timed sound onset on MacOS/X, and MS-Windows (ASIO requried)Low latency sound on Windows needs The standard sound subsystem of Microsoft Windows is not capable of reliably timed sound playback and

31、 capture with low latency If you want low latency sound with high timing precision on Microsoft Windows, youll need a sound card with native ASIO driver support and a special plugin for Psychtoolbox 采用ASIO(AudioStream Input Output)技术，可以减少系统对音频流信号的延迟，增强声卡硬件的处理能力。Delta 1010 声卡声卡Low latency,timed sound

32、 with PsychPortAudioSimple interface, allows for asynchronous sound output at ascheduled system time, e.g., time tvisualonset:1. Padevice = PsychPortAudio(Open, deviceid, mode, 2, 96000);2. Mysound = 0.9 * MakeBeep(1000, 0.1, 96000);3. PsychPortAudio(FillBuffer, Padevice, Mysound);4. PsychPortAudio(

33、Start, Padevice, 5, tvisualonset);5. Visonset = Screen(Flip, window, tvisualonset 0.004);6. .whatever.7. Audioonset = PsychPortAudio(Stop, Padevice);8. PsychPortAudio(Close ,Padevice); Auto-selects settings for low latency, but overrides possible. High precision, low latency with standard hardware o

34、n OS/X.PsychPortAudio-ASIO driverInitializePsychSound(1);freq = 96000;ISOI = 0.12; % 120 msStimdur=0.03; % auditory durationlatbias = (64 / freq); %hardware delay pahandle = PsychPortAudio(Open, ,2,freq);% Tell driver about hardwares inherent latencyprelat = PsychPortAudio(LatencyBias, pahandle, lat

35、bias);postlat = PsychPortAudio(LatencyBias, pahandle);%1: single sound; 2: two sounds; 3: multiple soundstone20 = 0.9 * MakeBeep(500,0.02,freq); %tone with 20 mstone(1,:) = 0.9 * MakeBeep(500, ISOI+Stimdur , freq); tone(2,1:length(tone20) = tone20; tone(2,end-length(tone20)+1:end) = tone20;for i=0:2

36、:6 range = 1+i*length(tone20): length(tone20)*(i+1); tone(3, range) = tone20;end PsychPortAudio(FillBuffer, pahandle, repmat(tone(1,:),2,1); PsychPortAudio(Start, pahandle, 1, 0); % PsychPortAudio(Stop, pahandle);使用前装载声使用前装载声卡驱动程序卡驱动程序List10.4Timingtry w=Screen(OpenWindow,0); WaitSecs(10); % 等待时间，以秒

37、为单位等待时间，以秒为单位 Screen(CloseAll);catch Screen(Close, w) rethrow(lasterror)endtry Screen(OpenWindow,0); t0= GetSecs; WaitSecs(5); t1= GetSecs; Screen(CloseAll ); t_elapsed = t1-t0; （反应时）（反应时）catch Screen(Close, w); rethrow(lasterror)endtry screenNum=0;res=1280 1024;clrdepth=32;wPtr,rect=Screen(OpenWind

38、ow,screenNum,0,0 0 res(1) res(2), clrdepth);black=BlackIndex(wPtr);white=WhiteIndex(wPtr);Screen(FillRect,wPtr,black);priorityLevel = MaxPriority(wPtr);Priority(priorityLevel);Screen(wPtr, Flip);Priority(0); Screen(CloseAll);catch Screen(Close, wPtr) rethrow(lasterror)EndPriority() - Switch Matlab p

39、rocess to realtime-scheduling mode.Priority() - Switch Matlab process to realtime-scheduling mode. 优先级优先级Basic commands for system control and timingT = GetSecs （获取时间值）（获取时间值） Query time with microsecond resolution. Uses highest resolution system clock for measurement of time.WaitSecs(duration) Wait

40、 for a specified amount of time duration. On MacOS/X: Accurate to 0.3 milliseconds on average. On MS-Windows: Accurate to KbDemoKbWaitdisp(Testing KbWait: hit any key. Just once.);startSecs = GetSecs;timeSecs = KbWait;keyIsDown, t, keyCode = KbCheck; str=KbName(keyCode), typed at time, num2str(timeS

41、ecs - startSecs), secondsdisp(str);Key Presses Collecting keypresses using pause and input Inputdlg promptParameters = Subject Name, Age, Gender (F or M?),Handedness (L or R);defaultParameters = , , R;Subinfo = inputdlg(promptParameters, Subject Info , 1, defaultParameters); Input resp_num= input (p

42、ress a number key); resp_char = input (press a number key);disp(Using input command)resp=x;while resp=a & resp=bresp=input(press a or b . , s);endrespBear in mind that the command window needs to be in front for the key-press to be available to Matlab, so input doesnt work well if you are using Psyc

43、htoolbox. In that case you will need to use GetChar, CharAvail, or KbCheck instead.% CharAvailWaitSecs(1)disp(Using CharAvail command)disp (Wait 2 sec to see if you pressed a key during that time);FlushEventsticwhile toc0 break; end end answer = find(buttons,1);MouseHideCursor;ShowCursor;SetMouse;Ge

44、tClicks List10.14 List10.15List10.16dioIn=digitalio(parallel,LPT1);dioOut=digitalio(parallel,LPT1);addline(dioOut,0:7,out);addline(dioIn,11:12,in); getvalue(dioIn); if sum(getvalue(dioIn)=1 response=0; elseif sum(getvalue(dioIn)=0 response=1; elseif sum(getvalue(dioIn)=2 response=2; end Keyboard/Mou

45、se Manipulations Adjustment with keyboard- Kingdom (1997, Perception)-see KeyboardMani.m Mller-Lyer illusion (Mller-Lyer, 1889) see mousemanipulation.mList10.17 Psychtoolbox实验设计 内容 实验设计步骤和编程概要 刺激的随机化 文件操作、结果保存 实验设计实例 GUI图形界面实验步骤 实验前关闭屏幕保护程序。 环境设置（图形界面） 指导语的呈现 输入被试个人资料（图形界面） 进行实验 记录结果和数据初步整理刺激的随机化 re

46、shape fullfact rand/seed randperm gentrials mod 被试内、被试间和重复次数应用举例 视觉搜索(Visual Search) 适应后效(Adaptation aftereffect) 协同运动(Coherent motion) 变化盲(Change Blindness) 快速序列视觉呈现(RSVP) 视听整合：audiovisual Ternus, what you see is what you hear, bouncing illusion1、视觉搜索Set size: 4 levels(4,8,12,16)Conjunction vs. fea

47、ture search;Target present vs. absent参考visualSearch, visualsearch2011文件夹，以及Ma et al. (2011)-Nature Neurosci程序框架和函数runVisualSearch.msettingsVisualSearch.mmakeSearchDisplay.m2、 适应后效-Adaptation to natural facial categories (Nature 2004：428, 558-561)程序框架-主函数和调用函数faceAdaptationMain.mgetParams.mdoFaceAdap

48、t.mdoFaceAdaptTrial.mopenFaceAdpt.mshowImage.mgetResponse.mSmallBlank.mfunction showImage(imName, wPtr, dur)img=imread(imName, JPG);textureIndex=Screen(MakeTexture, wPtr, double(img);Screen(DrawTexture, wPtr, textureIndex);Screen(wPtr, Flip);WaitSecs(dur)function response=getResponse(params)FlushEve

49、nts;sound(params.beep);response=GetChar;function doFaceAdapt(params, wPtr, dur, time)lastImage=0; for repeat=1:(time/dur)whichImage=ceil(rand*params.numAdaptImages);while whichImage=lastImage whichImage=ceil(rand*params.numAdaptImages); endimName=params.adaptstr,num2str(whichImage);showImage(imName,

50、 wPtr, dur);lastImage=whichImage;end完整的程序clear all; close all;cd(/Users/ariel/PTBTutorial/FaceExperiment/)addpath(pwd);addpath(morph);params = getParams;response = ;Screen(Preference, VisualDebuglevel, 3)wPtr rect = openFaceAdapt;doFaceAdapt (params, wPtr, params.preAdaptImageDuration, params.preAda

51、ptTime);for trialCounter=1:params.nTrials responsetrialCounter = doFaceAdaptTrial (wPtr, params);if trialCounterparams.nTrials doFaceAdapt(params, wPtr, params.reAdaptImageDuration, params.reAdaptTime);endend Screen(CloseAll);ShowCursorfunction params=getParamsparams.preAdaptTime=10; params.preAdaptImageDuration=2; params.numAdaptImages=10; params.numTestImages=49; %There are 49 morph imagesparams.testImageDuration=0.5; %Presentation of the test image (in secs)params.reAdaptTime=6; %Top up adaptationpara