diff --git a/.gitignore b/.gitignore index 869ff7bc3b8757071a6b6de32660f338706e4edc..51870579848044cf722a94abcf919dc5fd446170 100644 --- a/.gitignore +++ b/.gitignore @@ -31,3 +31,14 @@ Ex3/ex3_task3a.flt Ex3/ex3_task3b.flt Ex3/ex3_task3b.wav Ex3/log_file.txt + +Ex3-solution + +Ex4/out_win_dbg +Ex4/out_win_rls +Ex4/out_linux_dbg +Ex4/out_linux_rls + +Ex4/Ex4_task2_dbg.exe +Ex4/Ex4_task2.dot +Ex4/log_file.txt diff --git a/Ex4/.vscode/c_cpp_properties.json b/Ex4/.vscode/c_cpp_properties.json new file mode 100644 index 0000000000000000000000000000000000000000..11c576f694c99f7b30d0d5015bc0fede255966c2 --- /dev/null +++ b/Ex4/.vscode/c_cpp_properties.json @@ -0,0 +1,65 @@ +{ + "configurations": [ + { + "name": "Windows-Release", + "includePath": [ + "${workspaceFolder}/**", + "${workspaceFolder}/../DSPE_lib_minGW/MinGW-W64_8.1.0/include/**" + ], + "defines": [ + "WIN32", + "__DEBUG__=0" + ], + //"compilerPath": "gcc.exe", + "cStandard": "c11", + "cppStandard": "c++11", + "intelliSenseMode": "gcc-x64" + }, + { + "name": "Linux-Release", + "includePath": [ + "${workspaceFolder}/**", + "${workspaceFolder}/../DSPE_lib_minGW/MinGW-W64_8.1.0/include/**" + ], + "defines": [ + "WIN32", + "__DEBUG__=0" + ], + "compilerPath": "gcc.exe", + "cStandard": "c11", + "cppStandard": "c++11", + "intelliSenseMode": "gcc-x64" + }, + { + "name": "Windows-Debug", + "includePath": [ + "${workspaceFolder}/**", + "${workspaceFolder}/../DSPE_lib_minGW/MinGW-W64_8.1.0/include/**" + ], + "defines": [ + "WIN32", + "__DEBUG__=1" + ], + //"compilerPath": "gcc.exe", + "cStandard": "c11", + "cppStandard": "c++11", + "intelliSenseMode": "gcc-x64" + }, + { + "name": "Linux-Debug", + "includePath": [ + "${workspaceFolder}/**", + "${workspaceFolder}/../DSPE_lib_minGW/MinGW-W64_8.1.0/include/**" + ], + "defines": [ + "WIN32", + "__DEBUG__=1" + ], + "compilerPath": "gcc.exe", + "cStandard": "c11", + "cppStandard": "c++11", + "intelliSenseMode": "gcc-x64" + } +], + "version": 4 +} \ No newline at end of file diff --git a/Ex4/.vscode/launch.json b/Ex4/.vscode/launch.json new file mode 100644 index 0000000000000000000000000000000000000000..1a11b370e9af53a852d931e0428bae4b7332f93d --- /dev/null +++ b/Ex4/.vscode/launch.json @@ -0,0 +1,29 @@ +{ + // Use IntelliSense to learn about possible attributes. + // Hover to view descriptions of existing attributes. + // For more information, visit: https://go.microsoft.com/fwlink/?linkid=830387 + "version": "0.2.0", + "configurations": [ + { + "name": "(gdb) Ex3 â widnows debug run ", + "type": "cppdbg", + "request": "launch", + "program": "${workspaceFolder}/Ex4_task2_dbg.exe", + "args": [], + "stopAtEntry": false, + "cwd": "${workspaceFolder}", + "environment": [], + "externalConsole": false, + "MIMode": "gdb", +// "miDebuggerPath": "d:/CodeBlocks_20_03/MinGW/bin/gdb.exe", + "setupCommands": [ + { + "description": "WĹÄ cz formatowanie kodu dla gdb", + "text": "-enable-pretty-printing", + "ignoreFailures": true + } + ], + "preLaunchTask": "Build Ex4_task2.cpp" + } + ] +} \ No newline at end of file diff --git a/Ex4/.vscode/tasks.json b/Ex4/.vscode/tasks.json new file mode 100644 index 0000000000000000000000000000000000000000..90f99ae144e58cb1d5bf0706c48a73de2e25a398 --- /dev/null +++ b/Ex4/.vscode/tasks.json @@ -0,0 +1,33 @@ +{ + // See https://go.microsoft.com/fwlink/?LinkId=733558 + // for the documentation about the tasks.json format + "version": "2.0.0", + // https://code.visualstudio.com/docs/editor/tasks + // ${command:cpptools.activeConfigName} + "tasks": [ + { + "label": "Build Ex4_task2.cpp", + "type": "shell", + "command": "make build -f Makefile.main FILE=Ex4_task2 VS_CFG=${command:cpptools.activeConfigName}", + "group": { + "kind": "build", + "isDefault": true + }, + "presentation": { + "echo": true + }, + "problemMatcher": "$gcc" + }, + { + "label": "Clean Ex4_task2.cpp", + "type": "shell", + "command": "make clean -f Makefile.main VS_CFG=${command:cpptools.activeConfigName}", + "group": "build", + "presentation": { + "echo": true + }, + "problemMatcher": "$gcc" + } + + ] +} \ No newline at end of file diff --git a/Ex4/Ex4_task2.cbp b/Ex4/Ex4_task2.cbp new file mode 100644 index 0000000000000000000000000000000000000000..2ee317b6091f1944b25a700603ee6b20c105d34b --- /dev/null +++ b/Ex4/Ex4_task2.cbp @@ -0,0 +1,38 @@ +<?xml version="1.0" encoding="UTF-8" standalone="yes" ?> +<CodeBlocks_project_file> + <FileVersion major="1" minor="6" /> + <Project> + <Option title="Ex4_task2" /> + <Option pch_mode="2" /> + <Option compiler="gcc" /> + <Build> + <Target title="default"> + <Option output="Ex4_task2" prefix_auto="1" extension_auto="1" /> + <Option type="1" /> + <Option compiler="gcc" /> + </Target> + </Build> + <Compiler> + <Add option="-std=c++0x" /> + <Add option="-m32" /> + <Add option="-g" /> + <Add option="-DWIN32" /> + <Add directory="../DSPE_lib_minGW/MinGW-W64_8.1.0/include" /> + <Add directory="../DSPE_lib_minGW/MinGW-W64_8.1.0/include/dbg" /> + </Compiler> + <Linker> + <Add option="-static-libgcc" /> + <Add option="-m32" /> + <Add library="DSPE" /> + <Add library="winmm" /> + <Add directory="../DSPE_lib_minGW/MinGW-W64_8.1.0/dbg" /> + </Linker> + <Unit filename="Ex4_task2.cpp" /> + <Extensions> + <code_completion /> + <envvars /> + <debugger /> + <lib_finder disable_auto="1" /> + </Extensions> + </Project> +</CodeBlocks_project_file> diff --git a/Ex4/Ex4_task2.cpp b/Ex4/Ex4_task2.cpp new file mode 100644 index 0000000000000000000000000000000000000000..d47708a469bba8db3e77a921c898fd764f996909 --- /dev/null +++ b/Ex4/Ex4_task2.cpp @@ -0,0 +1,128 @@ +/*! Laboratory: Advanced signal processing of digital telecommunications + * (Zaawansowane przetwarzanie sygnaĹĂłw telekomunikacji cyfrowej) + * Ex. 4. task 2. + * Sampling rate conversion with fractional delay filters + * + * \author Marek Blok + * \date 2021.06.29 + */ +#include <DSP_lib.h> + +int main(int argn, char *args[]) +{ + /*************************************************************/ + // Log file setup + DSP::log.SetLogFileName("log_file.txt"); + DSP::log.SetLogState(DSP::e::LogState::file | DSP::e::LogState::console); + + DSP::log << DSP::lib_version_string() << endl << endl; + /*************************************************************/ + +/* + Idea: + - set of filters is read from a file - L filters designed in MATLAB + - tests with bandlimited or fullband signals - chirps + spectrograph + + 1. Load table of impulse responses of FSD filters + + 2. Input -> output buffer with callback + + DSP::u::OutputBuffer (unsigned int BufferSize_in, unsigned int NoOfInputs_in, DSPe_buffer_type cyclic, DSP_clock_ptr ParentClock, DSP_clock_ptr NotificationsClock, unsigned int NoOfOutputs_in, DSP::u::buffer_callback_ptr func_ptr, unsigned int CallbackIdentifier=0) + If NoOfOutputs_in > 0 NotificationsClock has also the meaning of OutputClock. + + 3. Synchronous clocks >> although itmight be better if the asynchronous clocks were used + +*/ + + + DSP::LoadCoef coef_info; + long int Fp1, Fp2; + unsigned int L, M; + unsigned int N_all; + DSP::Float_vector h_all; + + if (coef_info.Open("ex4_task2_h_FSD_all.coef", "matlab") == false) + { + return -1; + } + Fp1 = coef_info.Fp; + + // verify L and calculate M + DSP::Float_vector tmp; + coef_info.Load(tmp, 1); + L = (unsigned int)tmp[0]; M = (unsigned int)tmp[1]; + Fp2 = (Fp1*L)/M; + + // ********************************** + N_all = coef_info.GetSize(0); + coef_info.Load(h_all, 0); + + /*************************************************************/ + + DSP::Clock_ptr InputClock; + InputClock=DSP::Clock::CreateMasterClock(); + + + DSP::u::FileInput InputSignal(InputClock, "matlab/test_signal.wav", 1U, DSP::e::SampleType::ST_short, DSP::e::FileType::FT_wav); + int Fp1_tmp = InputSignal.GetSamplingRate(); + if (Fp1_tmp != Fp1) + { + DSP::log << DSP::e::LogMode::Error + << "Sampling rate of the input signal is different from sampling rate from fitler coefficients file (ex4_task2_h_FSD_all.coef)" << endl; + } + + DSP::u::SamplingRateConversion FSDresampler(InputClock, L, M, h_all); + + DSP::u::AudioOutput SoundOut(Fp2, 1, 16); + DSP::u::FileOutput FileOut_a("ex4_task2.wav", DSP::e::SampleType::ST_short, 1, DSP::e::FileType::FT_wav, Fp2); + DSP::u::FileOutput FileOut_b("ex4_task2.flt", DSP::e::SampleType::ST_short, 1, DSP::e::FileType::FT_flt, Fp2); + + DSP::log << "Fp1 = " << Fp1 << ", Fp2 = " << Fp2 << ", L = " << L << ", M = " << M << endl; + + /*************************************************************/ + // Connections definitions + InputSignal.Output("out") >> FSDresampler.Input("in"); + FSDresampler.Output("out") >> SoundOut.Input("in"); + FSDresampler.Output("out") >> FileOut_a.Input("in"); + FSDresampler.Output("out") >> FileOut_b.Input("in"); + + + ///////////////////////////////// + // check if there are signals + // connected to all inputs + DSP::Component::CheckInputsOfAllComponents(); + + // *********************************** // + DSP::Clock::SchemeToDOTfile(InputClock, "Ex4_task2.dot"); + + // *********************************** // + int SamplesInSegment = 512; + __int64 NoOfSamplesProcessed = 0; + // 10 seconds + #define MAX_SAMPLES_TO_PROCESS 10*Fp1 + while(NoOfSamplesProcessed < MAX_SAMPLES_TO_PROCESS) + { + + // ********************************************************** // + DSP::Clock::Execute(InputClock, SamplesInSegment); + // ********************************************************** // + + if (InputSignal.GetBytesRead() > 0) + { + NoOfSamplesProcessed = 0; // Play the whole file + } + else // Play 200ms more + { + if (NoOfSamplesProcessed < MAX_SAMPLES_TO_PROCESS - Fp1/5) + NoOfSamplesProcessed = MAX_SAMPLES_TO_PROCESS - Fp1/5; + } + + NoOfSamplesProcessed += SamplesInSegment; + // ********************************************************** // + } + + /*************************************************************/ + DSP::Clock::FreeClocks(); + + return 0; +} diff --git a/Ex4/Makefile b/Ex4/Makefile new file mode 100644 index 0000000000000000000000000000000000000000..c94c79f51791f3dffcb970c87be83384cef43774 --- /dev/null +++ b/Ex4/Makefile @@ -0,0 +1,85 @@ +# Run: make Release +# Run: make Debug +CC=g++ +# comflag = -m32 +# comflag = -m64 +comflag = $(COMFLAG) + +DSPElib_DIR = ../../_DSPE_lib_minGW_/$(DSPElib_SUBDIR) + +SRC_DIR = . +SRC_CPP_SUBDIR = . + +#DFLAGS = -DWIN32 -DDEVCPP + +# -D INCLUDE_DSPE_EXAMPLES # TODO: uĹźycie w ramach kompilacji Main.cpp w trybie DEBUG +ifeq ($(MODE),Release) + CFLAGS = $(comflag) -std=c++0x -O3 -Wall -c -fmessage-length=0 -fno-strict-aliasing + LINKER_FLAGS = $(comflag) -s -static-libgcc -static-libstdc++ $(MISC_LINKER_FLAGS) + INCLUDES := -I"$(DSPElib_DIR)/include" -I"$(DSPElib_DIR)/include/rls" + DSPElib_FULLDIR = $(DSPElib_DIR)/rls + EXE_FILENAME = $(CPP_FILENAME)_rls.exe +else + CFLAGS = $(comflag) -std=c++0x -O0 -g3 -Wall -c -fmessage-length=0 -W -Wshadow -Wconversion -fstrict-aliasing -fmax-errors=5 + LINKER_FLAGS = $(comflag) -static-libgcc -static-libstdc++ $(MISC_LINKER_FLAGS) + INCLUDES := -I"$(DSPElib_DIR)/include" -I"$(DSPElib_DIR)/include/dbg" + DSPElib_FULLDIR = $(DSPElib_DIR)/dbg + EXE_FILENAME = $(CPP_FILENAME)_dbg.exe +endif +# -U__STRICT_ANSI__ jest potrzebne do kompilacji debug_new.cpp, jezeli pominac ten plik to mozna rowniez wyrzucic te opcje +#CFLAGS_debug = $(comflag) -std=c++0x -O0 -g3 -Wall -c -fmessage-length=0 -W -Wshadow -Wco#nversion -fstrict-aliasing -U__STRICT_ANSI__ + +SOURCES_NAMES = +SOURCES_NAMES += $(CPP_FILENAME).cpp +SOURCES = $(addprefix $(SRC_CPP_SUBDIR)/,$(SOURCES_NAMES)) + +SOURCES_DBG = +# SOURCES_DBG += $(SRC_DIR)/Main.cpp + +# ################################################# # +# DEBUG +OBJECTS := $(SOURCES:%.cpp=$(OUT_DIR)/%.o) +DEPENDS := $(SOURCES:%.cpp=$(OUT_DIR)/%.d) + +# ################################################# # +-include $(DEPENDS) + +all: build + + +# ########################################################################################### # +# ########################################################################################### # +build: $(SRC_DIR)/$(EXE_FILENAME) + +$(SRC_DIR)/$(EXE_FILENAME): $(OBJECTS) + @echo $(EXE_FILENAME) + $(CC) -L$(DSPElib_FULLDIR) $(OBJECTS) -o"$(SRC_DIR)/$(EXE_FILENAME)" $(LINKER_FLAGS) -lDSPE $(LIBS) + +# ########################################################################################### # +# ########################################################################################### # +# Z podanej listy usuwany $(OUT_DIR_WIN_RLS)/ oraz '.o' zamieniamy na '.cpp' +$(OBJECTS): $(OUT_DIR)/%.o : %.cpp + @echo $(@D) $< $@ + + #mkdir -p $(OUT_DIR)/$(SRC_CPP_SUBDIR) + mkdir -p $(@D) + $(CC) $(DFLAGS) $(CFLAGS) $(INCLUDES) -MMD $< -o $@ + + +clean: + @echo MODE: $(MODE) + + @if [ -d "$(OUT_DIR)" ]; then \ + echo "cleaning $(OUT_DIR_DBG) ..."; \ + #find $(OUT_DIR)/ -name "*.o" -type f -delete; \ + rm -rf $(OUT_DIR)/*.d; \ + rm -rf $(OUT_DIR)/*.o; \ + rm -rf $(OUT_DIR); \ + echo "cleaned $(OUT_DIR)"; \ + fi + rm -rf "$(SRC_DIR)/$(EXE_FILENAME)"; \ + #rm -rf "$(SRC_DIR)/*.gif"; \ + #rm -rf "$(SRC_DIR)/*.dot"; \ + +.PHONY: all build clean + diff --git a/Ex4/Makefile.main b/Ex4/Makefile.main new file mode 100644 index 0000000000000000000000000000000000000000..cf6c368f484646f2b2a80c30863e62c3d1e7ec7f --- /dev/null +++ b/Ex4/Makefile.main @@ -0,0 +1,60 @@ +# (View > Command Palette) => "Convert Indentation to Tabs" + +ifeq ($(VS_CFG),Windows-Debug) + MAKEFILE = "Makefile" + MODE = Debug + COMFLAG = -m64 + + OUT_DIR = ./out_win_dbg + DSPElib_SUBDIR = MinGW-W64_8.1.0 + MISC_LINKER_FLAGS = -static + LIBS = -lwinmm -lws2_32 + DFLAGS = -DWIN32 -DDEVCPP +endif +ifeq ($(VS_CFG),Windows-Release) + MAKEFILE = "Makefile" + MODE = Release + COMFLAG = -m64 + + OUT_DIR = ./out_win_rls + DSPElib_SUBDIR = MinGW-W64_8.1.0 + MISC_LINKER_FLAGS = -static + LIBS = -lwinmm -lws2_32 + DFLAGS = -DWIN32 -DDEVCPP +endif +ifeq ($(VS_CFG),Linux-Debug) + MAKEFILE = "Makefile" + MODE = Debug + COMFLAG = + + OUT_DIR = ./out_linux_dbg + DSPElib_SUBDIR = $(shell gcc -dumpmachine)-gcc_$(shell gcc -dumpversion) + MISC_LINKER_FLAGS = + LIBS := -lasound + DFLAGS = +endif +ifeq ($(VS_CFG),Linux-Release) + MAKEFILE = "Makefile" + MODE = Release + COMFLAG = + + OUT_DIR = ./out_linux_rls + DSPElib_SUBDIR = $(shell gcc -dumpmachine)-gcc_$(shell gcc -dumpversion) + MISC_LINKER_FLAGS = + LIBS := -lasound + DFLAGS = +endif + + + +build: + @echo "Building $(VS_CFG)" + @echo $(VS_CFG): $(MODE) // $(MAKEFILE) + make build CPP_FILENAME=$(FILE) MODE=$(MODE) COMFLAG=$(COMFLAG) DFLAGS="$(DFLAGS)" LIBS="$(LIBS)" OUT_DIR=$(OUT_DIR) DSPElib_SUBDIR=$(DSPElib_SUBDIR) MISC_LINKER_FLAGS="$(MISC_LINKER_FLAGS)" -f $(MAKEFILE) + +clean: + @echo "Cleaning $(VS_CFG)" + @echo $(VS_CFG): $(MODE) // $(MAKEFILE) + make clean MODE=$(MODE) OUT_DIR=$(OUT_DIR) -f $(MAKEFILE) + + diff --git a/Ex4/matlab/FSDfilter.m b/Ex4/matlab/FSDfilter.m new file mode 100644 index 0000000000000000000000000000000000000000..7b941ff43f3d1064430e8f83da10de9e134b32d4 --- /dev/null +++ b/Ex4/matlab/FSDfilter.m @@ -0,0 +1,142 @@ +function [h, PE_FSD] = FSDfilter(method, N, delay, voff_typ) +% delay - total filter delay + +switch method + case 'LSE' + h=LSE(N,delay); + case 'LSE2' + h=LSE(N,delay,voff_typ); + case 'Lagr' + h=Lagr(N,delay); + case 'window' + h=window(N,delay,voff_typ); + case 'chebwin' + h=window(N,delay,3,voff_typ); + case 'offset' + h=offset(N,delay,voff_typ); + case 'optimal' + eval('h=getFSD(N,voff_typ,delay-(N-1)/2);','h=NaN;'); + + if isnan(h) + Krok=0; + Factor=4; + + while isnan(h) + Krok=Krok+1; + krok=Krok; voff=voff_typ; + while krok + voff=voff+(1/2-voff)/Factor; + krok=krok-1; + end; + eval('[h, PE_FSD, ekst]=getFSD(N,voff,delay-(N-1)/2);','h=NaN;'); + h=h.'; + end; + + while Krok>=0 + krok=Krok; voff=voff_typ; + while krok + voff=voff+(1/2-voff)/Factor; + krok=krok-1; + end; + ekst=ekst/ekst(end)*voff; + [h, PE_FSD, ekst]=getFSD(N,voff,delay-(N-1)/2, ekst); + h=h.'; + Krok=Krok-1; + end; + end; + + case 'optimal2' + % rasterless implementation + [h, PE_FSD] =getFSD_new(N,voff_typ,delay-(N-1)/2); +% eval('h=getFSD_new(N,voff_typ,delay-(N-1)/2);','h=NaN;'); + + if isnan(h) + Krok=0; + Factor=4; + + while isnan(h) + Krok=Krok+1; + krok=Krok; voff=voff_typ; + while krok + voff=voff+(1/2-voff)/Factor; + krok=krok-1; + end; + eval('[h, PE_FSD, ekst]=getFSD_new(N,voff,delay-(N-1)/2);','h=NaN;'); + h=h.'; + end; + + while Krok>=0 + krok=Krok; voff=voff_typ; + while krok + voff=voff+(1/2-voff)/Factor; + krok=krok-1; + end; + ekst=ekst/ekst(end)*voff; + [h, PE_FSD, ekst]=getFSD_new(N,voff,delay-(N-1)/2, ekst); + h=h.'; + Krok=Krok-1; + end; + end; +end; + +function h=Lagr(N,delay) + +for n=0:N-1 + m=0:N-1; + m(n+1)=[]; + h(n+1)=prod((delay-m)./(n-m)); +end; + +function h=LSE(N,delay, voff) + +n=0:N-1; +if nargin ==2, + h=sinc(n-delay); +elseif nargin==3 + alfa=2*voff; + for k=0:N-1 + l=0:N-1; + P(k+1,l+1)=alfa*sinc(alfa*(k-l)); + p(k+1,1)=alfa*sinc(alfa*(k-delay)); + end; + h=inv(P)*p; + h = h(:).'; +end; + +function h=window(N,delay,typ, Rp) +n=0:N-1; +hsinc=sinc(n-delay); + +if typ==3 + wsym=chebwin(N,Rp).'; +else + switch typ + case 2 + C=[0.43 -0.5 0.08]; %Blackman + case 1 + C=[0.54 -0.46 0]; %Hamming + case 0 + C=[0.5 -0.5 0]; %Hann + end; + wsym=C(1)+C(2)*cos(2*pi*n/(N-1))+C(3)*cos(4*pi*n/(N-1)); +end; + +h=hsinc.*wsym; + +function h=offset(N,delay,typ) +n=0:N-1; +hsinc=sinc(n-delay); +switch typ + case 2 + C=[0.43 -0.5 0.08]; %Blackman + case 1 + C=[0.54 -0.46 0]; %Hamming + case 0 + C=[0.5 -0.5 0]; %Hann +end; + +epsilon=delay-(N-1)/2; +%wsym=C(1)+C(2)*cos(2*pi*n/(N-1))+C(3)*cos(4*pi*n/(N-1)); +woffset=C(1)+C(2)*cos(2*pi*(n-epsilon)/(N-1))+C(3)*cos(4*pi*(n-epsilon)/(N-1)); +h=hsinc.*woffset; + diff --git a/Ex4/matlab/SPECgraf.m b/Ex4/matlab/SPECgraf.m new file mode 100644 index 0000000000000000000000000000000000000000..d892ec7febd085494a1abaf1f3c254b455018b45 --- /dev/null +++ b/Ex4/matlab/SPECgraf.m @@ -0,0 +1,1436 @@ +function varargout = SPECgraf(Akcja, param) +%cos(cumsum(n/1000)) + +% \fixed 2017.03.22 psd ==> pwelch +% \fixed 2012.03.05 finite ==> isfinite +% \fixed 2006.10.11 dealt with LineStyle/Marker warning +% \fixed 2005.11.20 dealt with "Warning: Log of zero." in ylim evaluation +% \fixed 2005.11.20 Placed all in single file +% \fixed 2005.11.20 Fixed problems with zoom out + +if nargin == 0 % LAUNCH GUI + +% fig = openfig(mfilename,'reuse'); return; %ZapiszFig(1,'specGUI.m') + if isempty(findobj(0, 'tag', 'Specgraf_DrawFig')) + + fig = Specgraf_DrawFig; + set(fig,'tag', 'Specgraf_DrawFig', 'Units', 'pixels'); + set(fig,'name', 'Spektrogram sygnałów niestacjonarnych (2017.03.22) dr inż. Marek Blok)',... + 'KeyPressFcn','1;'); + + SPECgraf('Init'); + SPECgraf('signal'); + set(fig,'Visible','on'); +% SPECgraf('per'); + else + SPECgraf('Exit'); + end + return; +end; + +if ~isstr(Akcja) % INVOKE NAMED SUBFUNCTION OR CALLBACK + disp('Something''s wrong'); + return; +end; + +% Generate a structure of handles to pass to callbacks, and store it. +fig=findobj(0, 'tag', 'Specgraf_DrawFig'); +hEditN=findobj(fig,'tag', 'N_edit'); +hEditM=findobj(fig,'tag', 'M_edit'); +hEditL=findobj(fig,'tag', 'L_edit'); +hEditK=findobj(fig,'tag', 'K_edit'); +hEditO=findobj(fig,'tag', 'O_edit'); +hEditW=findobj(fig,'tag', 'w_edit'); +hEditX=findobj(fig,'tag', 'x_edit'); +h_real=findobj(fig,'tag', 'real_checkbox'); +h_dB=findobj(fig,'tag', 'dB_checkbox'); +hEditNoise=findobj(fig,'tag', 'Noise_edit'); +hEditName=findobj(fig,'tag', 'Name_edit'); +hMenu=findobj(fig,'tag', 'choose_popupmenu'); +h_det=findobj(fig,'tag', 'detrend_checkbox'); +ha=get(fig, 'UserData'); +Ktory=get(hMenu,'Value'); +hEdit_dY=findobj(fig,'tag', 'dY_edit'); + +if strcmp(Akcja, 'Exit') + close(fig); + return; +elseif strcmp(Akcja, 'Init') + set(hEditX,'UserData','randn(1,L)'); + set(hEditX,'String','randn(1,L)'); + set(hEditX,'Max', 2); + + set(hEditW,'UserData','boxcar(M)'); + set(hEditW,'String','boxcar(M)'); + + set(hEditL,'UserData','1000'); + set(hEditL,'String','1000'); + set(hEditM,'UserData','100'); + set(hEditM,'String','100'); + set(hEditNoise,'UserData','-100'); + set(hEditNoise,'String','-100'); + +% set(hEditN,'UserData','1'); +% set(hEditN,'String','1'); + set(hEditO,'UserData','0'); + set(hEditO,'String','0'); + + set(hEditK,'UserData','256'); + set(hEditK,'String','256'); + + set(hEditName,'String','new'); + set(hMenu,'String','new'); + + set(h_real,'UserData',1); + set(h_real,'Value',1); + + set(h_dB,'UserData',0); + set(h_dB,'Value',0); + + set(h_det,'UserData',1); + set(h_det,'Value',1); + + ha(1)=findobj(fig,'tag', 'Signal_re_axes'); + ha(2)=findobj(fig,'tag', 'Signal_im_axes'); + ha(3)=findobj(fig,'tag', 'spec_axes'); + ha(4)=findobj(fig,'tag', 'per_axes'); +% for ind=[1 2 4], axes(ha(ind)); zoom on; end; +% axes(ha(3)); zoom off; + + set(fig, 'UserData', ha); + + set(hMenu,'UserData', [NaN, NaN, NaN, NaN, NaN, NaN, NaN, 1, -100, 5, -100, 5]); %handles & maximal values + set(hEdit_dY,'String','120'); + return; +elseif strcmp(Akcja, 'change_name') + pom=get(hMenu,'String'); + pom2=get(hEditName,'String'); + pom(Ktory,1:length(pom2)+1)=[pom2 0]; + set(hMenu,'String',pom); + set(hMenu, 'Value', Ktory); + return; +elseif strcmp(Akcja, 'new') + pom=get(hMenu,'String'); + Ktory=size(pom,1)+1; + pom(Ktory,1:4)=['new' 0]; + set(hMenu,'String',pom); + + + pom=get(hEditX,'UserData'); + pom(Ktory,1:11)=['randn(1,L)' 0]; + set(hEditX,'UserData',pom); + pom=get(hEditW,'UserData'); + pom(Ktory,1:10)=['boxcar(M)' 0]; + set(hEditW,'UserData',pom); + + pom=get(hEditL,'UserData'); + pom(Ktory,1:4)=['100' 0]; + set(hEditL,'UserData',pom); + pom=get(hEditM,'UserData'); + pom(Ktory,1:4)=['100' 0]; + set(hEditM,'UserData',pom); + pom=get(hEditNoise,'UserData'); + pom(Ktory,1:5)=['-100' 0]; + set(hEditNoise,'UserData',pom); + +% pom=get(hEditN,'UserData'); +% pom(Ktory,1:2)=['1' 0]; +% set(hEditN,'UserData',pom); + pom=get(hEditO,'UserData'); + pom(Ktory,1:2)=['0' 0]; + set(hEditO,'UserData',pom); + + pom=get(hEditK,'UserData'); + pom(Ktory,1:5)=['256' 0]; + set(hEditK,'UserData',pom); + + pom=get(hMenu,'UserData'); + pom(Ktory,1:size(pom,2))=ones(1,size(pom,2))*NaN; + set(hMenu,'UserData',pom); + + pom=get(h_real,'UserData'); + pom(Ktory,1)=1; + set(h_real,'Value', pom(Ktory,1)); + set(h_real,'UserData',pom); + + pom=get(h_det,'UserData'); + pom(Ktory,1)=1; + set(h_det,'Value', pom(Ktory,1)); + set(h_det,'UserData',pom); + + set(hMenu,'Value', Ktory); + + SPECgraf('Choose'); + SPECgraf('signal', Ktory); + return; +elseif strcmp(Akcja, 'delete') + pom=get(hMenu,'String'); + if size(pom,1)==1, + %SPECgraf('Reset'); + return; + end; + pom(Ktory,:)=[]; + set(hMenu,'String',pom); + + + pom=get(hEditX,'UserData'); + pom(Ktory,:)=[]; + set(hEditX,'UserData',pom); + pom=get(hEditW,'UserData'); + pom(Ktory,:)=[]; + set(hEditW,'UserData',pom); + + pom=get(hEditL,'UserData'); + pom(Ktory,:)=[]; + set(hEditL,'UserData',pom); + pom=get(hEditM,'UserData'); + pom(Ktory,:)=[]; + set(hEditM,'UserData',pom); + pom=get(hEditNoise,'UserData'); + pom(Ktory,:)=[]; + set(hEditNoise,'UserData',pom); + +% pom=get(hEditN,'UserData'); +% pom(Ktory,:)=[]; +% set(hEditN,'UserData',pom); + pom=get(hEditO,'UserData'); + pom(Ktory,:)=[]; + set(hEditO,'UserData',pom); + + pom=get(hEditK,'UserData'); + pom(Ktory,:)=[]; + set(hEditK,'UserData',pom); + + pom=get(hMenu,'UserData'); + for ind=1:3, + if isfinite(pom(Ktory,ind)) + delete(pom(Ktory,ind)); + end + end + pom(Ktory,:)=[]; + set(hMenu,'UserData',pom); + + pom=get(h_real,'UserData'); + pom(Ktory,:)=[]; + set(h_real,'UserData',pom); + + pom=get(h_det,'UserData'); + pom(Ktory,:)=[]; + set(h_det,'UserData',pom); + + set(hMenu,'Value', 1); + +% SPECgraf('signal', Ktory); + SPECgraf('Choose'); + return; + + +elseif strcmp(Akcja, 'Choose') + %selected new filter response + pom=get(hMenu,'String'); + ind=find(pom(Ktory,:)==0); + if isempty(ind) + ind=size(pom,2); + else + ind=ind(1)-1; + end + set(hEditName,'String',pom(Ktory,1:ind)); + + pom=get(hEditX,'UserData'); + set(hEditX, 'String', pom(Ktory,:)); + pom=get(hEditW,'UserData'); + set(hEditW, 'String', pom(Ktory,:)); + + pom=get(hEditL,'UserData'); + set(hEditL, 'String', pom(Ktory,:)); + pom=get(hEditM,'UserData'); + set(hEditM, 'String', pom(Ktory,:)); + pom=get(hEditNoise,'UserData'); + set(hEditNoise, 'String', pom(Ktory,:)); + +% pom=get(hEditN,'UserData'); +% set(hEditN, 'String', pom(Ktory,:)); + pom=get(hEditO,'UserData'); + set(hEditO, 'String', pom(Ktory,:)); + + pom=get(hEditK,'UserData'); + set(hEditK, 'String', pom(Ktory,:)); + + pom=get(h_real,'UserData'); + set(h_real,'Value', pom(Ktory,1)); + + pom=get(h_det,'UserData'); + set(h_det,'Value', pom(Ktory,1)); + +% SPECgraf('signal', Ktory); + return; + +elseif strcmp(Akcja, 'signal') + if nargin==2, + Ktory=param; + else + %save strings + tekst=[get(hEditL, 'String') 0]; + pom=get(hEditL,'UserData'); + pom(Ktory,1:max([1; size(pom,2); length(tekst)]))= ... + zeros(1, max([1; size(pom,2); length(tekst)])); + pom(Ktory,1:length(tekst))=tekst; + set(hEditL,'UserData', setstr(pom)); + + tekst=[get(hEditNoise, 'String') 0]; + pom=get(hEditNoise,'UserData'); + pom(Ktory,1:max([1; size(pom,2); length(tekst)]))= ... + zeros(1, max([1; size(pom,2); length(tekst)])); + pom(Ktory,1:length(tekst))=tekst; + set(hEditNoise,'UserData', setstr(pom)); + + tekst=get(hEditX, 'String').'; + tekst=[tekst(:); 0].'; + pom=get(hEditX,'UserData'); + pom(Ktory,1:max([1; size(pom,2); length(tekst)]))= ... + zeros(1, max([1; size(pom,2); length(tekst)])); + pom(Ktory,1:length(tekst))=tekst; + set(hEditX,'UserData', setstr(pom)); + + tekst=get(h_real, 'Value'); + pom=get(h_real,'UserData'); + pom(Ktory,1)=tekst; + set(h_real,'UserData', pom); + end; + + pom=get(hMenu,'UserData'); + hp_re=pom(:,1); + hp_im=pom(:,2); + hp_spec_t=pom(:,3); + hp_spec_t2=pom(:,4); + hp_spec=pom(:,5); + hp_per=pom(:,6); + hp_per2=pom(:,7); + max_x=pom(:,8); + min_spec=pom(:,9); + max_spec=pom(:,10); + min_per=pom(:,11); + max_per=pom(:,12); + + %generate signal + tekstL=get(hEditL,'UserData'); tekstL=tekstL(Ktory,:); + ind=find(tekstL==0); + if ~isempty(ind) + tekstL=tekstL(1:ind(1)-1); + end; + eval(['L=' tekstL ';'], 'L=1;') + + n=0:L-1; + tekstX=get(hEditX,'UserData'); tekstX=tekstX(Ktory,:); + ind=find(tekstX==0); + if ~isempty(ind) + tekstX=tekstX(1:ind(1)-1); + end; + eval(['x=' tekstX ';'], 'x=zeros(1,L);') + + Re=get(h_real,'UserData'); Re=Re(Ktory,1); + + tekstNoise=get(hEditNoise,'UserData'); tekstNoise=tekstNoise(Ktory,:); + ind=find(tekstNoise==0); + if ~isempty(ind) + tekstNoise=tekstNoise(1:ind(1)-1); + end; + eval(['Noise=' tekstNoise ';'], 'Noise=-300;') + + x=x(:); + if length(x)<L; + x(L)=0; + else + x=x(1:L); + end; + + N_lin=10.^(Noise/20); + if Re==1 + x=real(x)+N_lin*randn(size(x)); + else + x=x+N_lin*(randn(size(x))+j*randn(size(x)))/sqrt(2); + end; + max_x(Ktory)=max(abs([real(x); imag(x)])); + + %draw signal + + %real part of the signal + axes(ha(1)); + if isnan(hp_re(Ktory)) + hold on; + hp_re(Ktory)=plot(n,real(x), 'Color', 'b'); + hold off; + else + set(hp_re(Ktory),'Xdata', n, 'Ydata', real(x)); + end +% set(ha(1), 'Xlim', [-0.5, L-0.5], 'Ylim', [-1.1 1.1]*max(max_x)); +% eval('zoom reset', 'set(get(ha(1),''ZLabel''),''UserData'',[]);'); +% reset(get(ha(1),'ZLabel')); + + %imaginary part of the signal +% axes(ha(2)); + if isnan(hp_im(Ktory)) + hold on; + hp_im(Ktory)=plot(n,imag(x), 'Color', 'r'); + hold off; + else + set(hp_im(Ktory),'Xdata', n, 'Ydata', imag(x)); + end + set(ha(1), 'Xlim', [-0.5, L-0.5], 'Ylim', [-1.1 1.1]*max(max_x)); +% set(get(ha(2),'ZLabel'),'UserData',[]); +% reset(get(ha(2),'ZLabel')); +% eval('zoom reset', 'set(get(ha(1),''ZLabel''),''UserData'',[]);'); + if L>512 + set([hp_re, hp_im], 'Marker', '.', 'MarkerSize', 4); + else + set([hp_re, hp_im], 'LineStyle', '-'); + end; + eval('rmappdata(get(ha(1),''Zlabel''),''ZOOMAxesData'')','set(get(ha(1),''ZLabel''),''UserData'',[])'); + + set(hMenu,'UserData', [hp_re, hp_im, hp_spec_t, hp_spec_t2, hp_spec, hp_per, hp_per2, max_x, min_spec, max_spec, min_per, max_per]); + + %compute and draw periodogram + SPECgraf('spec', Ktory) +% SPECgraf zoom_on; + return; + +elseif strcmp(Akcja, 'spec') + if nargin==2, + Ktory=param; + else + %save strings + tekst=[get(hEditK, 'String') 0]; + pom=get(hEditK,'UserData'); + pom(Ktory,1:max([1; size(pom,2); length(tekst)]))= ... + zeros(1, max([1; size(pom,2); length(tekst)])); + pom(Ktory,1:length(tekst))=tekst; + set(hEditK,'UserData', setstr(pom)); + + tekst=[get(hEditM, 'String') 0]; + pom=get(hEditM,'UserData'); + pom(Ktory,1:max([1; size(pom,2); length(tekst)]))= ... + zeros(1, max([1; size(pom,2); length(tekst)])); + pom(Ktory,1:length(tekst))=tekst; + set(hEditM,'UserData', setstr(pom)); + +% tekst=[get(hEditN, 'String') 0]; +% pom=get(hEditN,'UserData'); +% pom(Ktory,1:max([1; size(pom,2); length(tekst)]))= ... +% zeros(1, max([1; size(pom,2); length(tekst)])); +% pom(Ktory,1:length(tekst))=tekst; +% set(hEditN,'UserData', setstr(pom)); + + tekst=[get(hEditO, 'String') 0]; + pom=get(hEditO,'UserData'); + pom(Ktory,1:max([1; size(pom,2); length(tekst)]))= ... + zeros(1, max([1; size(pom,2); length(tekst)])); + pom(Ktory,1:length(tekst))=tekst; + set(hEditO,'UserData', setstr(pom)); + + tekst=[get(hEditW, 'String') 0]; + pom=get(hEditW,'UserData'); + pom(Ktory,1:max([1; size(pom,2); length(tekst)]))= ... + zeros(1, max([1; size(pom,2); length(tekst)])); + pom(Ktory,1:length(tekst))=tekst; + set(hEditW,'UserData', setstr(pom)); + + tekst=get(h_det, 'Value'); + pom=get(h_det,'UserData'); + pom(Ktory,1)=tekst; + set(h_det,'UserData', pom); + end; + + pom=get(hMenu,'UserData'); + hp_re=pom(:,1); + hp_im=pom(:,2); + hp_spec_t=pom(:,3); + hp_spec_t2=pom(:,4); + hp_spec=pom(:,5); + hp_per=pom(:,6); + hp_per2=pom(:,7); + max_x=pom(:,8); + min_spec=pom(:,9); + max_spec=pom(:,10); + min_per=pom(:,11); + max_per=pom(:,12); + + %generate signal + tekstK=get(hEditK,'UserData'); tekstK=tekstK(Ktory,:); + ind=find(tekstK==0); + if ~isempty(ind) + tekstK=tekstK(1:ind(1)-1); + end; + eval(['K=' tekstK ';'], 'K=16;') + + tekstM=get(hEditM,'UserData'); tekstM=tekstM(Ktory,:); + ind=find(tekstM==0); + if ~isempty(ind) + tekstM=tekstM(1:ind(1)-1); + end; + eval(['M=' tekstM ';'], 'M=16;') + if M>K + M=K; + set(hEditM,'String', num2str(M)); + end; + +% tekstN=get(hEditN,'UserData'); tekstN=tekstN(Ktory,:); +% ind=find(tekstN==0); +% if ~isempty(ind) +% tekstN=tekstN(1:ind(1)-1); +% end; +% eval(['N=' tekstN ';'], 'N=16;') + + tekstO=get(hEditO,'UserData'); tekstO=tekstO(Ktory,:); + ind=find(tekstO==0); + if ~isempty(ind) + tekstO=tekstO(1:ind(1)-1); + end; +% eval(['O=' tekstO ';'], 'O=0;') + O=eval(tekstO, '0'); % \Fixed 2005.11.03 + O=round(O/100*M); % \Fixed 2005.11.03 nakładkowanie podawane w procentach !!! + + + tekstW=get(hEditW,'UserData'); tekstW=tekstW(Ktory,:); + ind=find(tekstW==0); + if ~isempty(ind) + tekstW=tekstW(1:ind(1)-1); + end; + eval(['w=' tekstW ';'], 'w=ones(1,M);') + + w=w(:); + if length(w)<M; + w(M)=0; + else + w=w(1:M); + end; + + x=get(hp_re(Ktory), 'Ydata')+j*get(hp_im(Ktory), 'Ydata'); + +% O=floor(O/100*M); % \Fixed 2005.11.03 + if O>=M + O=M-1; + set(hEditO,'String', num2str(O/M*100)); + end; + N = fix((length(x)-O)/(M-O)); + set(hEditN, 'String', sprintf('%i', N)); + +% det_=get(h_det,'UserData'); det_=det_(Ktory,1); +% if det_==0, +% det='none'; +% else +% det_='linear'; +% end; + + [Spec, f, t]=specgram(x, K, 1, w, O); + Spec=(abs(Spec).^2)/norm(w)^2; +% [Per, f2]=psd(x, K, 1, w, O); + [Per, f2]=pwelch(x, w, O, K, 1); + Spec=abs(Spec); + + Re=get(h_real,'UserData'); + if ~Re(Ktory,1) + f=fftshift(f); + f(1:floor(K/2))=f(1:floor(K/2))-1; + f2=fftshift(f2); + f2(1:floor(K/2))=f2(1:floor(K/2))-1; + Spec=[Spec(ceil(K/2):K,:); Spec(1:floor(K/2),:)]; + Per=fftshift(Per); + end + + min_spec(Ktory)=min([min(Spec(isfinite(Spec))); 0]); + max_spec(Ktory)=max([max(Spec(isfinite(Spec))); 0.001]); + min_per(Ktory)=min([min(Per(isfinite(Per))); 0]); + max_per(Ktory)=max([max(Per(isfinite(Per))); 0.001]); + if get(h_dB, 'UserData') %dB + Spec=10*log10(Spec); + Per=10*log10(Per); + end + set(get(ha(2),'Ylabel'),'UserData', f); + set(get(ha(3),'Ylabel'),'UserData', Spec); + set(get(ha(4),'Ylabel'),'UserData', Per); + + %draw signal + axes(ha(2)); + if length(t)>1 + t2=t+(t(2)-t(1))/2; + else + t2=M/2; + end; + + if isnan(hp_spec_t(Ktory)) + hold on; +% hp_spec_t(Ktory)=plot(t2,max(Spec), 'Color', 'k'); + pomoc=abs(get(hp_re(Ktory),'Ydata')+j*get(hp_im(Ktory),'Ydata')); + hp_spec_t(Ktory)=plot(get(hp_re(Ktory),'Xdata'), pomoc, 'Color', 'k'); + if length(pomoc)>512 + set(hp_spec_t(Ktory), 'Marker', '.', 'MarkerSize', 4); + else + set(hp_spec_t(Ktory), 'LineStyle', '-'); + end; + hold off; + else +% set(hp_spec_t(Ktory),'Xdata', t2, 'Ydata', max(Spec), 'Color', 'k'); + set(hp_spec_t(Ktory),'Xdata', get(hp_re(Ktory),'Xdata'),... + 'Ydata', abs(get(hp_re(Ktory),'Ydata')+j*get(hp_im(Ktory),'Ydata')), 'Color', 'k'); + end + if length(t)==1, + set(ha(2), 'Xlim', [t(1)-0.5 t(1)+0.5], 'Ylim', [-1.1 1.1]*max(max_x)); + else + set(ha(2), 'Xlim', [t(1) t(length(t))+(t(2)-t(1))], 'Ylim', [-1.1 1.1]*max(max_x)); + end; +% set(get(ha(2),'ZLabel'),'UserData',[]); +% reset(get(ha(2),'ZLabel')); + eval('rmappdata(get(ha(2),''Zlabel''),''ZOOMAxesData'')','set(get(ha(2),''ZLabel''),''UserData'',[])'); + + %spektrogram + axes(ha(3)); + if isnan(hp_spec(Ktory)) + hold on; + hp_spec(Ktory)=image(t2, f, Spec); + colormap(hot); + hold off; + else + set(hp_spec(Ktory),'Xdata', t2, 'Ydata', f, 'Cdata', Spec); + end + if length(t)==1, + tlim_=[t(1)-0.5 t(1)+0.5]; + else + tlim_=[t(1) t(length(t))+t(2)-t(1)]; + end; + if all(Re) + set(ha(3), 'Ylim', [0 0.5], 'Xlim', tlim_); + else + set(ha(3), 'Ylim', [-0.5 0.5], 'Xlim', tlim_); + end +% set(get(ha(3),'ZLabel'),'UserData',[]); +% reset(get(ha(3),'ZLabel')); + eval('zoom reset', 'set(get(ha(3),''ZLabel''),''UserData'',[]);'); + + axes(ha(4)); + if isnan(hp_per(Ktory)) + hold on; + hp_per(Ktory)=plot(Per, f2, 'Color', 'k'); + hold off; + else + set(hp_per(Ktory),'Ydata', f2, 'Xdata', Per, 'Color', 'k'); + end + set(ha(4), 'Xdir', 'reverse', 'YTick', []); %, 'Xlim', [t(1) t(length(t))], 'Ylim', [-1.1 1.1]*max(max_x)); + if all(Re) + set(ha(4), 'Ylim', [0 0.5]); + else + set(ha(4), 'Ylim', [-0.5 0.5]); + end + eval('rmappdata(get(ha(4),''Zlabel''),''ZOOMAxesData'')','set(get(ha(4),''ZLabel''),''UserData'',[])'); + + set(hMenu,'UserData', [hp_re, hp_im, hp_spec_t, hp_spec_t2, hp_spec, hp_per, hp_per2, max_x, min_spec, max_spec, min_per, max_per]); + + % delete + if 1 + pom=get(hMenu,'UserData'); + hp_=pom(:,[4,7]); + hp_=hp_(isfinite(hp_)); + if length(hp_)>0 + delete(hp_) + pom(:,4)=NaN; + pom(:,7)=NaN; + set(hMenu,'UserData',pom); + end; + end + + SPECgraf('spec_ylim'); + SPECgraf zoom_on; + SPECgraf zoom_off; + return; + +elseif strcmp(Akcja, 'dB') + pom=get(h_dB, 'UserData'); + if pom~=get(h_dB, 'Value') + pom=get(hMenu,'UserData'); + hp_spec=pom(:,3); + for ind=1:length(hp_spec) + sygn=get(get(ha(3),'Ylabel'),'UserData'); + per=get(get(ha(4),'Ylabel'),'UserData'); + if get(h_dB, 'Value') %dB + sygn=10*log10(sygn); + per=10*log10(per); + else %lin + sygn=10.^(sygn/10); + per=10.^(per/10); + end; + set(get(ha(3),'Ylabel'),'UserData', sygn); + set(get(ha(4),'Ylabel'),'UserData', per); +% set(hp_spec(ind),'Cdata', sygn); + end; + set(h_dB, 'UserData', get(h_dB, 'Value')); + + hp_=pom(:,[4 7]); + hp_=hp_(isfinite(hp_)); + if length(hp_)>0 + delete(hp_) + pom(:,4)=NaN; + pom(:,7)=NaN; + set(hMenu,'UserData',pom); + end; + SPECgraf('spec_ylim'); + end + return + +elseif strcmp(Akcja, 'spec_ylim') + pom=get(hMenu,'UserData'); + hp_re=pom(:,1); + hp_im=pom(:,2); + hp_spec_t=pom(:,3); + hp_spec=pom(:,5); + hp_per=pom(:,6); + min_spec=pom(:,9); + max_spec=pom(:,10); + min_per=pom(:,11); + max_per=pom(:,12); + if get(h_dB, 'UserData') %dB + tekst=get(hEdit_dY,'String'); + eval(['dY=' tekst ';'], 'dY=120;'); + if dY<=0, dY=10; end; + params_=[min(min_spec) max(max_spec)]; + ind_params = find(abs(params_) <= eps); + if length(ind_params) > 0, + params_(ind_params) = NaN*ones(size(ind_params)); + end + ylim_=10*log10(params_); + if ~isfinite(ylim_(2)) + ylim_(2)=0; + end + ylim_(1)=ylim_(2)-dY; + params_=[min(min_per) max(max_per)]; + ind_params = find(abs(params_) <= eps); + if length(ind_params) > 0, + params_(ind_params) = NaN*ones(size(ind_params)); + end + ylim_per=10*log10(params_); + if ~isfinite(ylim_per(2)) + ylim_per(2)=0; + end + ylim_per(1)=ylim_per(2)-dY; + else + ylim_=[0 max(max_spec)]; + ylim_per=[0 max(max_per)]; + end + ylim_per(2)=max([ylim_per(2) ylim_(2)]); + f=get(get(ha(2),'Ylabel'),'UserData'); + Spec=get(get(ha(3),'Ylabel'),'UserData'); + Per=get(get(ha(4),'Ylabel'),'UserData'); + set(hp_spec_t(Ktory),'Ydata', abs(get(hp_re(Ktory),'Ydata')+j*get(hp_im(Ktory),'Ydata'))); + set(ha(2),'Ylim', ylim_+(ylim_(2)-ylim_(1))*[-0.1 0.1]); + Spec=64*(Spec-ylim_(1))/(ylim_(2)-ylim_(1)); + + set(hp_spec(Ktory),'Cdata', Spec); + set(hp_per(Ktory),'Xdata', Per, 'Ydata', f); + set(ha(4),'Xlim', ylim_per); + + if get(h_dB, 'UserData') %dB + set(hp_spec_t(Ktory),'Ydata', 20*log10(abs(get(hp_re(Ktory),'Ydata')+j*get(hp_im(Ktory),'Ydata')))); + else + set(hp_spec_t(Ktory),'Ydata', abs(get(hp_re(Ktory),'Ydata')+j*get(hp_im(Ktory),'Ydata'))); + end + + SPECgraf zoom_on; + eval('zoom reset', 'set(get(ha(3),''ZLabel''),''UserData'',[]);'); + SPECgraf zoom_off; + return +elseif strcmp(Akcja, 'zoom_on') + zoom on; +% pom=get(fig, 'WindowButtonDownFcn'); +% set(fig, 'WindowButtonDownFcn', 'SPECgraf zoom'); +% set(get(ha(3),'Xlabel'), 'UserData', pom); + return; +elseif strcmp(Akcja, 'zoom_off') + if get(findobj(fig,'tag','checkbox_zoom'),'Value') == 0, +% pom = get(get(ha(3),'Xlabel'), 'UserData'); +% set(fig, 'WindowButtonDownFcn', pom); + zoom off; + set(fig, 'WindowButtonDownFcn', 'SPECgraf zoom'); + set(get(ha(3),'Xlabel'), 'UserData', '1;'); + end + return; +elseif strcmp(Akcja, 'zoom_spec') + if get(findobj(fig,'tag','checkbox_zoom'),'Value') ~= 0, + Specgraf zoom_on; + else + Specgraf zoom_off; + end +elseif strcmp(Akcja, 'zoom') +% if strcmp(get(fig,'SelectionType'),'normal') | (gca~=ha(3)) + pause(0); + if (get(gco,'Parent')~=ha(3)) | get(findobj(fig,'tag','checkbox_zoom'),'Value') + eval(get(get(ha(3),'Xlabel'), 'UserData')); + elseif get(gco,'Parent')==ha(3) + pom=get(ha(3), 'CurrentPoint'); + f_=pom(1,2); t_=pom(1,1); + + pom_=get(hMenu,'UserData'); + hp_spec_t2=pom_(:,4); + hp_spec=pom_(:,5); + hp_per2=pom_(:,7); + f=get(hp_spec(Ktory), 'Ydata'); + ind_f=find(abs(f-f_)==min(abs(f-f_))); ind_f=ind_f(1); + t=get(hp_spec(Ktory), 'Xdata'); +% if length(t)>1, +% t=t+(t(2)-t(1))/2; +% end; + ind_t=find(abs(t-t_)==min(abs(t-t_))); + ind_t=ind_t(1); + set(findobj(fig,'tag', 'text_t_f'),... + 'ForegroundColor', 'r', 'String', sprintf('n=%i, f=%6.3f', t(ind_t),f(ind_f))); + + Spec=get(get(ha(3),'Ylabel'),'UserData'); + axes(ha(4)); + if length(Spec(:,ind_t))>length(f) + ind=find(f==0); + Spec(ind(1),:)=[]; + end; + if isnan(hp_per2(Ktory)) + hold on; + hp_per2(Ktory)=plot(Spec(:,ind_t),f, 'Color', 'r', 'LineWidth', 2); + hold off; + else + set(hp_per2(Ktory),'Xdata', Spec(:,ind_t), 'Ydata', f); + end +% pom=get(ha(4), 'Xlim'); +% pom(2)=max([pom(2) max(Spec(:,ind_t))]); +% set(ha(4),'Xlim', pom); + + axes(ha(2)); + if isnan(hp_spec_t2(Ktory)) + hold on; + hp_spec_t2(Ktory)=plot(t,Spec(ind_f,:), 'Color', 'r', 'LineWidth', 2); + hold off; + SPECgraf zoom_on; + SPECgraf zoom_off; + else + set(hp_spec_t2(Ktory),'Xdata', t, 'Ydata', Spec(ind_f,:)); + end + + pom_(Ktory,7)=hp_per2; + pom_(Ktory,4)=hp_spec_t2; + set(hMenu,'UserData',pom_); + end; + return; +end; + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +function fig_h=Specgraf_DrawFig +fig_h=figure; +set(fig_h, ... + 'Units', 'Normalized',... + 'Position', [0.1125 0.1188 0.7813 0.7448],... + 'Color', [1.0000 1.0000 1.0000],... + 'Name', 'Untitled',... + 'NumberTitle', 'off',... + 'Menu', 'none',... + 'Tag', 'figure1'... + ); +h=axes; +set(h, ... + 'Units', 'Normalized',... + 'Position', [0.0330 0.0881 0.2150 0.4867],... + 'Color', [1.0000 1.0000 1.0000],... + 'XColor', [0.0000 0.0000 0.0000],... + 'YColor', [0.0000 0.0000 0.0000],... + 'FontSize', 10,... + 'Box', 'on',... + 'Tag', 'per_axes'... + ); +h=axes; +set(h, ... + 'Units', 'Normalized',... + 'Position', [0.2900 0.6350 0.6990 0.1538],... + 'Color', [1.0000 1.0000 1.0000],... + 'XColor', [0.0000 0.0000 0.0000],... + 'YColor', [0.0000 0.0000 0.0000],... + 'FontSize', 10,... + 'Box', 'on',... + 'Tag', 'Signal_im_axes'... + ); +h=axes; +set(h, ... + 'Units', 'Normalized',... + 'Position', [0.2880 0.0867 0.7030 0.4867],... + 'Color', [1.0000 1.0000 1.0000],... + 'XColor', [0.0000 0.0000 0.0000],... + 'YColor', [0.0000 0.0000 0.0000],... + 'FontSize', 10,... + 'Box', 'on',... + 'Tag', 'spec_axes'... + ); +h=axes; +set(h, ... + 'Units', 'Normalized',... + 'Position', [0.2920 0.8266 0.6980 0.1580],... + 'Color', [1.0000 1.0000 1.0000],... + 'XColor', [0.0000 0.0000 0.0000],... + 'YColor', [0.0000 0.0000 0.0000],... + 'FontSize', 10,... + 'Box', 'on',... + 'Tag', 'Signal_re_axes'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'frame',... + 'Units', 'Normalized',... + 'Position', [0.0050 0.7972 0.2400 0.1930],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', '',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'frame2'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'frame',... + 'Units', 'Normalized',... + 'Position', [0.0050 0.5916 0.2400 0.2028],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', '',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'frame4'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.1190 0.6364 0.0280 0.0350],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', '[%]',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text23'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'checkbox',... + 'Units', 'Normalized',... + 'Position', [0.8720 0.0126 0.0660 0.0196],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'zoom',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'Specgraf zoom_spec',... + 'Tag', 'checkbox_zoom'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'pushbutton',... + 'Units', 'Normalized',... + 'Position', [0.0140 0.0098 0.0660 0.0378],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'Exit',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf Exit',... + 'Tag', 'pushbutton5'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.4640 0.0028 0.3540 0.0350],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', '',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text_t_f'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'edit',... + 'Units', 'Normalized',... + 'Position', [0.1750 0.0126 0.0750 0.0364],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'Edit Text',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf(''spec_ylim'')',... + 'Tag', 'dY_edit'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.1110 0.0154 0.0620 0.0308],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'dY [dB] =',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text21'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'checkbox',... + 'Units', 'Normalized',... + 'Position', [0.9270 0.5748 0.0500 0.0322],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'dB',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf(''dB'')',... + 'Tag', 'dB_checkbox'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'pushbutton',... + 'Units', 'Normalized',... + 'Position', [0.0660 0.4909 0.0600 0.0420],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'Delete',... + 'Value', 0,... + 'Visible', 'off',... + 'Callback', 'SPECgraf(''delete'')',... + 'Tag', 'pushbutton4'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'pushbutton',... + 'Units', 'Normalized',... + 'Position', [0.0120 0.4909 0.0480 0.0420],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'New',... + 'Value', 0,... + 'Visible', 'off',... + 'Callback', 'SPECgraf(''new'')',... + 'Tag', 'pushbutton3'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.1580 0.6769 0.0750 0.0364],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'Edit Text',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf(''spec'')',... + 'Tag', 'N_edit'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.1280 0.6825 0.0300 0.0266],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'N =',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text18'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'edit',... + 'Units', 'Normalized',... + 'Position', [0.0390 0.8014 0.0720 0.0378],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'Edit Text',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf(''signal'')',... + 'Tag', 'L_edit'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.0090 0.8028 0.0290 0.0266],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'L =',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text17'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'checkbox',... + 'Units', 'Normalized',... + 'Position', [0.1290 0.7301 0.1100 0.0252],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'detrend',... + 'Value', 0,... + 'Visible', 'off',... + 'Callback', 'SPECgraf(''spec'')',... + 'Tag', 'detrend_checkbox'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'edit',... + 'Units', 'Normalized',... + 'Position', [0.0400 0.6378 0.0750 0.0364],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'Edit Text',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf(''spec'')',... + 'Tag', 'O_edit'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.0110 0.6378 0.0280 0.0294],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'O =',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text16'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'edit',... + 'Units', 'Normalized',... + 'Position', [0.0520 0.6000 0.1870 0.0336],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'Edit Text',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf(''spec'')',... + 'Tag', 'w_edit'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.0130 0.5986 0.0400 0.0294],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'w[n] =',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text15'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'edit',... + 'Units', 'Normalized',... + 'Position', [0.1670 0.8042 0.0720 0.0364],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'Edit Text',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf(''signal'')',... + 'Tag', 'Noise_edit'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.1190 0.8042 0.0470 0.0266],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'Noise=',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text14'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'edit',... + 'Units', 'Normalized',... + 'Position', [0.0410 0.6769 0.0750 0.0364],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'Edit Text',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf(''spec'')',... + 'Tag', 'M_edit'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.0110 0.6825 0.0300 0.0266],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'M =',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text13'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'edit',... + 'Units', 'Normalized',... + 'Position', [0.1340 0.4909 0.1030 0.0392],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'Edit Text',... + 'Value', 0,... + 'Visible', 'off',... + 'Callback', 'SPECgraf(''change_name'')',... + 'Tag', 'Name_edit'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'popupmenu',... + 'Units', 'Normalized',... + 'Position', [0.0120 0.5413 0.2260 0.0294],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'Popup Menu',... + 'Value', 1,... + 'Visible', 'off',... + 'Callback', 'SPECgraf(''Choose'')',... + 'Tag', 'choose_popupmenu'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'checkbox',... + 'Units', 'Normalized',... + 'Position', [0.1160 0.9580 0.1110 0.0266],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'real signal',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf(''signal'')',... + 'Tag', 'real_checkbox'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.0100 0.7622 0.1470 0.0266],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'Spectrograf parameters',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text10'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'edit',... + 'Units', 'Normalized',... + 'Position', [0.0400 0.7259 0.0750 0.0364],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'Edit Text',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf(''spec'')',... + 'Tag', 'K_edit'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.0140 0.7329 0.0240 0.0252],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'K =',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text9'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'edit',... + 'Units', 'Normalized',... + 'Position', [0.0390 0.8462 0.2010 0.1077],... + 'BackGroundColor', [1.0000 1.0000 1.0000],... + 'String', 'Edit Text',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', 'SPECgraf(''signal'')',... + 'Tag', 'x_edit'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.0130 0.9245 0.0240 0.0252],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'x =',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text4'... + ); +h=uicontrol; +x=version; +if (x(1)>='5'), + set(h, ... + 'FontSize', 10); +end; +set(h, ... + 'Style', 'text',... + 'Units', 'Normalized',... + 'Position', [0.0120 0.9636 0.1010 0.0238],... + 'BackGroundColor', [0.9255 0.9137 0.8471],... + 'String', 'Testing signal',... + 'Value', 0,... + 'Visible', 'on',... + 'Callback', '',... + 'Tag', 'text2'... + ); diff --git a/Ex4/matlab/ex4_task1.m b/Ex4/matlab/ex4_task1.m new file mode 100644 index 0000000000000000000000000000000000000000..a38fbcab999c6afa180ab4a7651b18a82205f309 --- /dev/null +++ b/Ex4/matlab/ex4_task1.m @@ -0,0 +1,44 @@ +function ex4_task1 + + +N_FSD = 13; +tau_N = (N_FSD-1)/2; + +Fp1 = 8000; Fp2 = 10000; + +M = Fp1 / gcd(Fp1,Fp2) +L = Fp2 / gcd(Fp1,Fp2) + +eps_all = -linspace(-0.5+0.5/L, 0.5-0.5/L, L); + +figure(1) +ha1 = subplot(2,1,1); +h_all = zeros(1, N_FSD*L); +n_all = 0:length(h_all)-1; n_all = n_all - (length(n_all)-1)/2; +for ind = 1:L, + epsilon = eps_all(ind); + h = FSDfilter('Lagr', N_FSD, tau_N+epsilon); +% h = FSDfilter('optimal2', N_FSD, tau_N+epsilon, 0.4); + h_all(ind:L:end) = h; + set(stem(ha1, n_all, h_all), 'Marker', 'none'); + pause(0.1) + drawnow +end + +K = 16*8192; +F = linspace(0, L*Fp1, K); + +H_all = abs(fft(h_all/L, K)); +figure(11) +subplot(2,1,1) +plot(h_all); +subplot(2,1,2) +plot(F, 20*log10(H_all)); +set(gca, 'Xlim', [0, L*Fp1/2], 'Ylim', [-100, 3]) + + +dane.h{1} = h_all; +dane.h{2} = [L, M]; +dane.Fp = Fp1; +save_filter_coef('ex4_task2_h_FSD_all', dane, 1); + diff --git a/Ex4/matlab/ex4_task1_test.m b/Ex4/matlab/ex4_task1_test.m new file mode 100644 index 0000000000000000000000000000000000000000..ec655b84ea590aee6879e51ec24ed92da3219df2 --- /dev/null +++ b/Ex4/matlab/ex4_task1_test.m @@ -0,0 +1,39 @@ +function ex4_task1_test + + +N_FSD = 63; +tau_N = (N_FSD-1)/2; + +% Fp1 = 8000; Fp2 = 10000; +Fp1 = 16000; Fp2 = 11025; + +M = Fp1 / gcd(Fp1,Fp2) +L = Fp2 / gcd(Fp1,Fp2) + +eps_all = -linspace(-0.5+0.5/L, 0.5-0.5/L, L); + +h_all = zeros(1, N_FSD*L); +for ind = 1:L, + epsilon = eps_all(ind); + h = FSDfilter('optimal2', N_FSD, tau_N+epsilon, 0.45); +% h = FSDfilter('optimal2', N_FSD, tau_N+epsilon, 0.4); + h_all(ind:L:end) = h; +end + +K = 16*8192; +F = linspace(0, L*Fp1, K); + +H_all = abs(fft(h_all/L, K)); +figure(11) +subplot(2,1,1) +plot(h_all); +subplot(2,1,2) +plot(F, 20*log10(H_all)); +set(gca, 'Xlim', [0, L*Fp1/2], 'Ylim', [-100, 3]) + + +dane.h{1} = h_all; +dane.h{2} = [L, M]; +dane.Fp = Fp1; +save_filter_coef('ex4_task2_h_FSD_all', dane, 1); + diff --git a/Ex4/matlab/ex4_task2_h_FSD_all.coef b/Ex4/matlab/ex4_task2_h_FSD_all.coef new file mode 100644 index 0000000000000000000000000000000000000000..b12e93ef7292943bc38b57f86abbe339bd50a6ad Binary files /dev/null and b/Ex4/matlab/ex4_task2_h_FSD_all.coef differ diff --git a/Ex4/matlab/getFSD_new.m b/Ex4/matlab/getFSD_new.m new file mode 100644 index 0000000000000000000000000000000000000000..ab7d7fefff185c6a67c53fe85c2d50f13fd19d7c --- /dev/null +++ b/Ex4/matlab/getFSD_new.m @@ -0,0 +1,468 @@ +% %poszukiwanie filtru optymalnego o dowolnej długoci +%w oparciu o układ równań rzeczywistych +% w oparciu o phase rotated real error +% +%implementacja w oparciu o poszukiwanie ekstremów +%pomiędzy dwoma kolejnymi punktami ekstremalnymi o tym samym znaku +% +% Wariant bezrastrowy w użyciem interpolacji parabolicznej +% +%getFSD_new(N,voff,d); +% N - długoć odpowiedzi impulsowej +% voff - częstotliwoć unormowana końca pasma dobrej aproksymacji +% d - opónienie netto + +% 2011.03.01 +function [h, TPEopt, ekst, A, B]=getFSD_new(N, voff, d, K, ekst, force_end); + +if nargin < 6, + force_end = 0; +end + +if nargin==0, + N=4; + voff=0.35; + d=0.25; + K=2048; %DFT length +elseif nargin==1, + voff=0.35; + d=0.25; + K=2048; %DFT length +elseif nargin==2 + d=0.25; + K=2048; %DFT length +elseif nargin==3 + K=2048; %DFT length +end; + +% df_factor = 10; df_factor_scaling = 1.05; +df_factor = 20; df_factor_scaling = 1.0; +df_factor = 4; + +% % K = K * 8 * 16; +% K_ = 8192*2*64; + +% N_LPF(ind)=ceil((-20*log10(sqrt(dp*ds))-13)/(14.6*(vs-vp))+1); +% N_FSD_est=ceil(-1/8*(PE/(4*(1/2-voff))+3)); +PE_est=-(N*8+3)*(4*(1/2-voff)); +if (PE_est < -140) + warning('Peek error estimation based on impulse response length shows that N is so large that it might result in fatal design errors'); +elseif (PE_est < -200) + warning('Peek error estimation based on impulse response length shows that N is so large that it probably will result in fatal design errors'); +elseif (PE_est < -275) + warning('Peek error estimation based on impulse response length shows that N is so large that it must result in fatal design errors'); + pause +end + + +% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% bez kwantyzacji osi częstotliwoci +M_in=round(voff*K)+1; % voff=(M_in-1)/K; + +L=(N-1)/2; % bulk delay +H_L = Hideal(L, voff, K); H_L=conj(H_L); +Hid = Hideal(L+d, voff, K); + +TPEopt=-300; +if nargin<5 +% if rem(N,2)==0 +% v=linspace(0,voff,N/2+1); +% else +% v=linspace(voff/N,voff,(N+1)/2); +% end; +% alfa=0.95; +% alfa=0.975; + alfa=0.999; + if (voff<=0.47) + alfa=0.99; + end +% alfa=0.996; +% alfa=0.999; + if rem(N,2)==0 +% '0' + x=[0:(N/2-1)]; + v=[0 cumsum(alfa.^x)]; + v=v/v(end)*voff; + else + x=[0:((N-1)/2)]; + v=cumsum(alfa.^x); + v(1)=v(1)/2; + v=v/v(end)*voff; + end; + ekst=(v*K); +% ekst=round(v*K); +end + +nn=1:N; A=zeros(N+1,N+1); B=zeros(N+1,1); +if (rem(N,2)==0) %parzyste + N_2=N/2; + A(1:N_2,N+1)=zeros(N_2,1); + A(((N_2+1):N),N+1)=((-1).^(1:N_2).'); + A(N+1,(1:(N+1)))=ones(1,N+1); + B(N+1)=1; +else + N_2=(N+1)/2; + A(1:N_2,N+1)=-((-1).^(1:N_2).'); + A(((N_2+1):(N+1)),N+1)=0; +end + +% ile = 0; +koniec = 1; +while 1, +% ile = ile + 1 + %construct equations + if (rem(N,2)==0) %parzyste + v=ekst(2:length(ekst))/K; %N-parzyste + for k=1:N_2, + A(k,nn)=sin(2*pi*(nn-1-L)*v(k)); + B(k)=sin(2*pi*v(k)*d); + A((N_2+k),nn)=cos(2*pi*(nn-1-L)*v(k)); + B((N_2+k))=cos(2*pi*v(k)*d); + end + else + v=ekst(1:length(ekst))/K; %N-nieparzyste + for k=1:N_2, + A(k,nn)=sin(2*pi*(nn-1-L)*v(k)); + B(k)=sin(2*pi*v(k)*d); + A((N_2+k),nn)=cos(2*pi*(nn-1-L)*v(k)); + B((N_2+k))=cos(2*pi*v(k)*d); + end + end + + if any(~isfinite(B)), + disp('B isn''t finite'); + error + end + if any(~isfinite(A)), + disp('A isn''t finite'); + error + end + x=A\B; % x=inv(A)*B; + h=x(1:N).'; + if any(~isfinite(h)), + disp('h isn''t finite'); + plot(20*log10(abs(Error))); + error + end + + if any(~isfinite(h)), + disp('Impulse response isn''t finite'); + error + end + + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % Error calculation + H = fft(h,K); + H = H(1:M_in); + H_L = H_L(1:M_in); + Hid = Hid(1:M_in); + + if (M_in-1)/K < voff, + n = 0:length(h)-1; + H_off = sum(h.*exp(-i*2*pi*voff*n)); + Hid_off = exp(-i*2*pi*voff*(L+d)); + H_L_off = exp(+i*2*pi*voff*L); + H = [H, H_off]; + Hid = [Hid, Hid_off]; + H_L = [H_L, H_L_off]; + end + + Error=H-Hid; + Error=Error.*H_L; + + n = 0:length(h)-1; + H_off = sum(h.*exp(i*2*pi*voff*(n-L).*n)); + Hid_off = exp(-i*2*pi*voff*d); + + % wyznaczenie współczynnika rotacji błędu + % => korygujemy błšd tak, żeby na częstotliwoci + Error_v_off = H_off - Hid_off; +% if abs(Error(M_in)) > 0, + if abs(Error_v_off) > 0, + Err_rot = conj(Error_v_off)/abs(Error_v_off); +% Err_rot = conj(Error(M_in))/abs(Error(M_in)); + Error=Error.*Err_rot; + else + Err_rot = 1; + end + if any(~isfinite(Error)), + disp('Normalized error isn''t finite'); + error + end + + if force_end == 1; + TPEopt = 20*log10(abs(x(N+1))); + TPEopt(2) = 20*log10(max(abs(Error))); + return; + end +% plot(20*log10(abs(Error))); + + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % poszukiwanie położenia nowych ekstremów +% ekst=Ekstr(real(Error),ekst, N); +% ekst=Ekstr_new(h, ekst, Err_rot, d, N); + ekst_old = ekst; + ekst=Ekstr_new(h, ekst/K, d, N, K, df_factor)*K; + df_factor = round(df_factor * df_factor_scaling); + if df_factor > 200, + df_factor = 200; + end + + if any(~isfinite(ekst)), + disp('ekst isn''t finite'); + error + end + + % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + % wyznaczenie peek error + %TPE=max(20*log10(abs(Error(ekst+1)))); + %TPE_H=max(20*log10(abs(Error))); + TPE=20*log10(abs(x(N+1))); + if TPEopt<TPE + TPEopt=TPE; + else +% ekst + ekst_old - ekst; + koniec = koniec - 1; + if koniec == 0, + break; + end + end +end + +return + + +%======================================================== +function Hid = Hideal(delay, F, K) +% delay - opónienie ułamkowe projektowanego filtru +% F - unormowana częstotliwoć górna pasma akuratnej aproksymacji +% K - liczba pršżków transformaty +% Hid - charakterystyka idealna (K-punktowa) + +v=(0:(K-1))/K; v(K/2:K)=v(K/2:K)-1; +Hid = ones(1,K).*exp(-i*2*pi*v*delay); + +%======================================================== +function x=Ekstr(Err, old, N) +% jeżeli filtr parzystej długoci to zakładamy, że punkty ekstremalne w 0 i M_in (end) +% nie ulegajš zmianie +% jeżeli filtr o nieparzystej długoci to w zerze napewno nie ma punktu ekstremalnego +% oraz punkt ekstremalny w M_in (end) nie ulega zmianie + +delta=max(abs(Err(old+1))); + +%wyjciowy zestaw punktów ekstremalnych +temp=old+1; %parzyste +if rem(N,2)==1 + temp=[1 temp]; %nieparzyste +end + +min_max=0; +for ind=length(temp):-1:3 + pom=Err(temp(ind-2):temp(ind)); + if min_max==0 + [wart, ind_new]=find(pom==min(pom)); + else + [wart, ind_new]=find(pom==max(pom)); + end + if abs(wart)>delta + temp(ind-1)=temp(ind-2)+ind_new-1; + end + min_max=1-min_max; +end + +x=temp-1; +if rem(N,2)==1 + x=x(2:length(x)); +end + +function f_new = Ekstr_new(h, f_old, d, N, K, df_factor); +% jeżeli filtr parzystej długoci to zakładamy, że punkty ekstremalne w 0 i M_in (end) +% nie ulegajš zmianie +% jeżeli filtr o nieparzystej długoci to w zerze napewno nie ma punktu ekstremalnego +% oraz punkt ekstremalny w M_in (end) nie ulega zmianie +L = (N-1)/2; +n = 0:N-1; n=n(:); + +f_old_ = f_old; +foff = f_old(end); +H_off = sum(h(:).*exp(-i*2*pi*foff*(n-L))); +Hid_off = exp(-i*2*pi*foff*d); +Error_f_off = H_off - Hid_off; +if abs(Error_f_off) > 0, + Err_rot = conj(Error_f_off)/abs(Error_f_off); +else + Err_rot = 1; +end + +for ind = 1:length(f_old), + H_off(ind) = sum(h(:).*exp(j*2*pi*f_old(ind)*(n-L).*n)); +end +Hid_old = exp(-i*2*pi*f_old*d); + +Err_old = real((H_off - Hid_old)*Err_rot); + +% delta=max(abs(Err(old+1))); +delta=max(abs(Err_old)); + + +do_draw = 0; +if do_draw == 1, + % % f = linspace(0,0.5, 2048); + f = linspace(0,0.5, K); % 2011.02.28 + + E_ = 0; + H_ = 0; + for n_ = n.', + E_ = E_ + Err_rot*(-j)*2*pi* h(n_+1).*(n_-(L+d)).*exp(-j*2*pi*f*(n_-(L+d))); + H_ = H_ + Err_rot*h(n_+1).*exp(-j*2*pi*f*(n_-L)); + end + + H_id = Err_rot*exp(-j*2*pi*f*d); + + figure(5) + subplot(3,1,1) + plot(f, real(E_)); + hold on + plot(f, imag(E_), 'r'); + hold off + set(gca, 'Xtick', f_old_, 'Xgrid', 'on'); + subplot(3,1,2) + plot(f, abs(H_-H_id), 'k'); + hold on + plot(f, real(H_-H_id), 'b'); + plot(f, imag(H_-H_id), 'r'); + hold off + + subplot(3,1,3) + plot(f, 20*log10(abs(H_-H_id)), 'k'); +end + +if rem(N, 2) == 1, + f_old = [0, f_old]; % 0 must always be included +end +% find extremas +% leave point at v_off (maximum) +f_new = f_old(end); + +polarization = -1; +while length(f_old) > 2, + % find extremum between f_start and f_end + f_end = f_old(end); + f_start = f_old(end-2); + + fo = f_old(end-1); % old value + + df = min([(f_end-fo), (fo-f_start)]) / df_factor; + + fo_new = find_extremum(h, d, Err_rot, [fo-df, fo, fo+df], polarization, do_draw); + + f_new = [fo_new, f_new]; + f_old(end) = []; + + polarization = -polarization; +end +if rem(N, 2) == 0, + f_new = [f_old(1), f_new]; % include 0 +end + +% size(f_old_) +% size(f_new) + +if do_draw == 1, + subplot(3,1,2) + set(gca, 'Xtick', f_new, 'Xgrid', 'on'); + subplot(3,1,3) + set(gca, 'Xtick', f_new, 'Xgrid', 'on'); + + % f_new + % + pause +end + +function fo_new = find_extremum(h, d, Err_rot, fs, polarization, do_draw) + +N = length(h); +n = 0:N-1; n = n(:); +L = (N-1)/2; + +f = fs; +% f = [fs(1), (fs(1)+fs(2))/2, fs(2), (fs(2)+fs(3))/2, fs(3)]; +% f = linspace(fs(1), fs(3), 11); + +% E_ = 0; +H_ = 0; +for n_ = n.', +% E_ = E_ + Err_rot*(-j)*2*pi* h(n_+1).*(n_-(L+d)).*exp(-j*2*pi*f*(n_-(L+d))); + H_ = H_ + Err_rot*h(n_+1).*exp(-j*2*pi*f*(n_-L)); +end +H_id = Err_rot*exp(-j*2*pi*f*d); + +E_ = real(H_-H_id) * polarization; + +%f(x) = a*x^2+b*x+c +%delta = b^2 - 4*a*c +%f_ekstremum = -b/(2*a) + +% % y = polyval(a, x) +% y1 = a*x1^2 + b*x1 + c +% y2 = a*x2^2 + b*x2 + c +% y3 = a*x3^2 + b*x3 + c + +a = polyfit(f, E_, 2); + +fo_new = -a(2)/(2*a(1)); + +% +% if (fo_new < fs(1)) +% fo_new = fs(1); +% end +% if (fo_new > fs(3)) +% fo_new = fs(3); +% end + +if (fo_new < fs(1)) || (fo_new > fs(3)) + ind = find(max(E_) == E_); + fo_new = f(ind(1)); +else + H_new = 0; + for n_ = n.', + H_new = H_new + Err_rot*h(n_+1).*exp(-j*2*pi*fo_new*(n_-L)); + end + H_id_new = Err_rot*exp(-j*2*pi*fo_new*d); + E_new = real(H_new-H_id_new) * polarization; + + E_ = [E_, E_new]; + f = [f, fo_new]; + + ind = find(max(E_) == E_); + if (fo_new ~= f(ind(1))), + fo_new = f(ind(1)); + end + +end + +% % if fo_new < fs(1) +% % fo_new = fs(1); +% % end +% % if fo_new > fs(3) +% % fo_new = fs(3); +% % end +% % % fo_new +% % % fo_new; + +% do_draw = 1; +if do_draw == 1, + f_ = linspace(fs(1), fs(3), 51); + y = polarization * polyval(a, f_); + subplot(3,1,2) + hold on + plot(f_, y, 'm'); + hold off + max_y = max(abs(y)); + if max_y == 0, max_y = 1, end; + set(gca, 'Ylim', [-1.15, 1.15]*max_y); + % % pause + % +end diff --git a/Ex4/matlab/readaudiofile.m b/Ex4/matlab/readaudiofile.m new file mode 100644 index 0000000000000000000000000000000000000000..825a9a890d520cb59e3dda822abed693cc2fc262 --- /dev/null +++ b/Ex4/matlab/readaudiofile.m @@ -0,0 +1,215 @@ +function [x, Fs] = readaudiofile(filename, param) +% [x, Fs] = readaudiofile(filename, param) +% +% returns vector x of the size SIZE=[samples channels]. +% eg. x = x(:,1) + j*x(:,2); +% +% special uses: +% param == 'size': returns SIZE in x +% param == 'cplx': for stereo files returns x as a complex vector instead of matrix +% +% supported file types: +% *.flt +% *.wav +% *.tape +% last modification: 2021.03.29 +% Author: Marek Blok + +return_cplx = 0; +if nargin == 1, + param = inf; +else + if strcmp(param, 'cplx') == 1, + return_cplx = 1; + param = inf; + end +end + +ind = find(filename == '.'); +if length(ind) == 0, + file_type = 'w'; % *.wav +else + ind = ind(end); + + temp = filename(ind+1:end); + + file_type = 'u'; % unknown format + if strcmp(temp, 'flt') == 1, + file_type = 'f'; % two channel floting point + elseif strcmp(temp, 'wav') == 1, + file_type = 'w'; % *.wav + elseif strcmp(temp, 'tape') == 1, + file_type = 't'; % *.tape + end +end + +switch file_type, + case 'w', + if strcmp(param, 'size') == 1, + if exist('audioread','file') == 0 + x = wavread(filename, 'size'); + % siz = [samples channels]. + else + info = audioinfo(filename); + x = [info.TotalSamples, info.NumChannels]; + Fs = info.SampleRate; + end + else + if isfinite(param) + if exist('audioread','file') == 0 + [x, Fs] = wavread(filename, param); + else + if length(param) == 1 + param = [1, param]; + end + [x, Fs] = audioread(filename, param); + end + else + if exist('audioread','file') == 0 + [x, Fs] = wavread(filename); + else + [x, Fs] = audioread(filename); + end + end + end + + case 't' + plik = fopen(filename, 'rb'); + if plik == -1, + error('File does not exist !!!'); + end + + header.size = fread(plik, 1, 'uint32', 0) + 4; + + header.fname = char(fread(plik, 256, 'char', 0).'); + header.cfg_fname = char(fread(plik, 256, 'char', 0).'); + header.sw_rev = fread(plik, 1, 'uint32', 0); + header.hw_rev = fread(plik, 1, 'uint32', 0); + header.file_ = fread(plik, 1, 'uint32', 0); + header.tape_type = fread(plik, 1, 'uint32', 0); + header.start_time = fread(plik, 1, 'int32', 0); % time_t + header.end_time = fread(plik, 1, 'int32', 0); % time_t + + header.total_samples = fread(plik, 1, 'uint32', 0); + file_length = header.total_samples * 4 + header.size + header.current_sample = fread(plik, 1, 'uint32', 0); + header.loop_start = fread(plik, 1, 'int64', 0); + header.loop_end = fread(plik, 1, 'int64', 0); + header.loop = fread(plik, 1, 'int32', 0); + header.group_size_32 = fread(plik, 1, 'uint32', 0); + header.block_size = fread(plik, 1, 'uint32', 0); + header.block_count = fread(plik, 1, 'uint32', 0); + header.fifo_size = fread(plik, 1, 'uint32', 0); + + + header.comment = char(fread(plik, 256, 'char', 0).'); + header.tmp = char(fread(plik, 20, 'char', 0).'); % time_t + header.status = fread(plik, 1, 'uint32', 0); + header.timestamps = fread(plik, 1, 'int32', 0); + header.freq = fread(plik, 1, 'float', 0); + header.cplx_datarate = fread(plik, 1, 'float', 0); + +% ftell(plik) + header.reserved = fread(plik, 128, 'uint32', 0); +% header.reserved.' + + header + ftell(plik) + + header.sample_type = 2; + header.ch_no = 2; + header.Fs = NaN; + + sample_type = 'int16'; + sample_size = 2; + + header_size = header.size; + if strcmp(param, 'size') == 1, + fseek(plik, 0, 'eof'); + size = (ftell(plik) - header_size) / sample_size / header.ch_no; % sizeof(float) *2 + x = size; + else + fseek(plik, header_size, 'bof'); + + len = param(1); + if length(param) > 1, + fseek(plik, sample_size*header.ch_no*(param(1)-1), 'cof'); + len = param(2) - param(1) + 1; + end + +% x = fread(plik, [header.ch_no, len], sample_type); + x = fread(plik, [header.ch_no, len], sample_type, 0); + x = x.'; + end + fclose(plik); + + case 'f' + plik = fopen(filename, 'rb'); + if plik == -1, + error('File does not exist !!!'); + end + + % 3 B - wersja pliku (w tym typ próbek) + % 1 B - liczba kanałów + % 4 B - szybkoć próbkowania + header_size = 8; + header.ver = fread(plik, 1, 'uint8'); + header.sample_type = fread(plik, 1, 'uint16'); + header.ch_no = fread(plik, 1, 'uint8'); + header.Fs = fread(plik, 1, 'uint32'); + + Fs = header.Fs; + + switch (header.ver), + case 0, + switch (header.sample_type), + case 0, + sample_type = 'float'; + sample_size = 4; + case 1, + sample_type = 'uchar'; + sample_size = 1; + case 2, + sample_type = 'short'; + sample_size = 2; + case 3, + sample_type = 'int'; + sample_size = 4; + otherwise + error('Unsupported *.flt sample type !!!'); + end + otherwise + error('Unsupported *.flt file version !!!'); + end + + + if strcmp(param, 'size') == 1, + fseek(plik, 0, 'eof'); + size = (ftell(plik) - header_size) / sample_size / header.ch_no; % sizeof(float) *2 + x = size; + else + len = param(1); + status = 0; + if length(param) > 1, + status = fseek(plik, sample_size*header.ch_no*(param(1)-1), 'cof'); + len = param(2) - param(1) + 1; + end + + if (status == -1) + x = []; + else + x = fread(plik, [header.ch_no, len], sample_type); + x = x.'; + end + end + + fclose(plik); + otherwise + error('Unsupported file format !!!'); +end + +if return_cplx == 1, + if length(x(1,:)) == 2, + x = x(:,1) + j*x(:,2); + end +end \ No newline at end of file diff --git a/Ex4/matlab/save_filter_coef.m b/Ex4/matlab/save_filter_coef.m new file mode 100644 index 0000000000000000000000000000000000000000..8cf1f072079de9ee5dab70dc9682ae0ae00ef1df --- /dev/null +++ b/Ex4/matlab/save_filter_coef.m @@ -0,0 +1,226 @@ +function save_filter_coef(filename, coefficients, file_version) +% save_filter_coef(filename, coefficients [, file_version]) +% +% FIR filter coefficients +% coefficients.h +% - can be cell vector to store multiple FIR impulse responces +% - if any of the impulse responses is complex then all responses will +% be stored as complex +% IIR filter coefficients +% coefficients.a +% coefficients.b +% +% For file_version >= 1 (default == 0), +% coefficients.Fp - sampling frequency +% +% This file is a part of Digital Signal Processing Engine +% +% File format (*.coef) - this is open format, for general use +% * (not only for storing coefficients) +% +% File format (*.h) - this is C++ header file for DSPElib with hardcoded *.coef equivalent +% +% \author Marek Blok +% \date 2020.02.22 + +if nargin == 2, + file_version = 0; +end + +ind = find(filename == '.'); +if length(ind) == 0, + filename = [filename, '.coef']; +end + +% detect mode based on file extention +[PATHSTR,NAME,EXT] = fileparts(filename); +switch EXT + case '.h' + mode = 'h-file'; + plik = fopen(filename, 'wt'); + + case '.coef' + mode = 'coefs-file'; + plik = fopen(filename, 'wb'); + + otherwise + mode = 'unknown-file'; + error(mode); +end + +switch mode + case 'h-file' + fprintf(plik, '// save_filter_coef output (hard coded coefficients)\n'); + fprintf(plik, '#pragma once\n'); + fprintf(plik, '#include <DSP_lib.h>\n'); + fprintf(plik, '\n'); + + if isfield(coefficients, 'h'), + h = coefficients.h; + + for ind_h = 0:length(h)-1 + fprintf(plik, 'DSP_float_vector Farrow_coefs_row_%i = {', ind_h); + h_tmp = h{ind_h+1}; + for ind_p = 0:length(h_tmp)-1, + if ind_p > 0 + fprintf(plik, ', '); + end + fprintf(plik, '%.15ff', h_tmp(ind_p+1)); + end + %0.1f, 0.2f, 0.1f + fprintf(plik, '};\n'); + end + fprintf(plik, '\n'); + + %vector <DSP_float_vector> Farrow_coefs = { coefs_0, coefs_1, coefs_2 }; + fprintf(plik, 'vector <DSP_float_vector> Farrow_coefs = {'); + for ind_h = 0:length(h)-1 + if ind_h > 0 + fprintf(plik, ', '); + end + if (ind_h < length(h_tmp)) & (rem(ind_h, 5) == 0) + fprintf(plik, '\n '); + end + fprintf(plik, 'Farrow_coefs_row_%i', ind_h); + end + fprintf(plik, '\n };\n'); + + fprintf(plik, '\n'); + fprintf(plik, 'const unsigned int p_Farrow = %i; // should be equal to Farrow_coefs[0].size()-1\n', length(h{1})-1); + fprintf(plik, 'const unsigned int N_FSD = %i; // should be equal to Farrow_coefs.size()\n', length(h)); + end + + if isfield(coefficients, 'b'), + fclose all; + error('unsupported: coefficients.b'); + end + if isfield(coefficients, 'a'), + fclose all; + error('unsupported: coefficients.a'); + end + + case 'coefs-file' + % * - (uchar) 1B - file version number + fwrite(plik, file_version, 'uchar'); + + switch file_version, + case 0, + case 1, + % * - (uint) 4B - Sampling frequency + if isfield(coefficients, 'Fp'), + fwrite(plik, coefficients.Fp, 'uint32'); + else + fclose(plik); + error('Input data does not contain Fp'); + end + otherwise, + fclose(plik); + error('This version of coefficients file is unsupported'); + end + + if isfield(coefficients, 'h'), + isFIR = 1; + if iscell(coefficients.h) + resp_no = length(coefficients.h); + % if resp_no = 1; + % coefficients.h = coefficients.h{1}; + % end + else + resp_no = 1; + end + if resp_no == 1, + isComplex = any(imag(coefficients.h(:))); + else + isComplex = false; + for ind_resp = 1:resp_no, + isComplex = isComplex | any(imag(coefficients.h{ind_resp}(:))); + end + end + else + isFIR = 0; + isComplex = any(imag(coefficients.a(:))) | any(imag(coefficients.b(:))); + end + + % * - data - coefficients data (depends on fle version) + % * . + % * Data segment format: + % * -# (version: 0x00) + % * - (uchar) 1B - number of sample dimensions + % * 1 - real, 2 - complex, ... + if isComplex, + fwrite(plik, 2, 'uchar'); % complex + else + fwrite(plik, 1, 'uchar'); % real + end + % * - (uchar) 1B - sample component type + % * - DSP_FT_float (=1) : C++ float (32bit floating point) + % * - DSP_FT_short (=2) : C++ short (16bit signed integer) + % * - DSP_FT_uchar (=3) : C++ unsigned char (8bit unsigned integer with bias (0x80)) + % * - DSP_FT_double (=7) : C++ double (64bit floating point) + % * - DSP_FT_long_double (=8) : C++ long double (80bit floating point) + fwrite(plik, 1, 'uchar'); + + % * - (uchar) 1B - number of vectors + % * - 1 - FIR filter coefficients (one vector) + % * - 2 - IIR filter coefficients (two vectors) + % * - (x number of vectors) + % * - (ushort) 2B - number of samples in vector + % * - (x number of samples) + % * - (x number of sample dimensions) + % * - (sample componet type) xB - sample component + % * e.g. real, imag part + if isFIR, + fwrite(plik, resp_no, 'uchar'); + + if iscell(coefficients.h) + for ind_resp = 1:resp_no, + N_FIR = length(coefficients.h{ind_resp}); + fwrite(plik, N_FIR, 'uint16'); + if isComplex, + dane(1:2:2*N_FIR) = real(coefficients.h{ind_resp}); + dane(2:2:2*N_FIR) = imag(coefficients.h{ind_resp}); + fwrite(plik, dane, 'float'); + else + fwrite(plik, real(coefficients.h{ind_resp}), 'float'); + end + end + else + N_FIR = length(coefficients.h); + fwrite(plik, N_FIR, 'uint16'); + if isComplex, + dane(1:2:2*N_FIR) = real(coefficients.h); + dane(2:2:2*N_FIR) = imag(coefficients.h); + fwrite(plik, dane, 'float'); + else + fwrite(plik, real(coefficients.h), 'float'); + end + end + + else + fwrite(plik, 2, 'uchar'); + + N_a = length(coefficients.a); + fwrite(plik, N_a, 'uint16'); + if isComplex, + dane(1:2:2*N_a) = real(coefficients.a); + dane(2:2:2*N_a) = imag(coefficients.a); + fwrite(plik, dane, 'float'); + else + fwrite(plik, real(coefficients.a), 'float'); + end + + + N_b = length(coefficients.b); + fwrite(plik, N_b, 'uint16'); + if isComplex, + dane(1:2:2*N_b) = real(coefficients.b); + dane(2:2:2*N_b) = imag(coefficients.b); + fwrite(plik, dane, 'float'); + else + fwrite(plik, real(coefficients.b), 'float'); + end + end + +end +fclose(plik); + diff --git a/Ex4/matlab/test1_11025.wav b/Ex4/matlab/test1_11025.wav new file mode 100644 index 0000000000000000000000000000000000000000..a47f814892bd68c7eb58e32d6c1502b956d20c9f Binary files /dev/null and b/Ex4/matlab/test1_11025.wav differ diff --git a/Ex4/matlab/test2_11025.wav b/Ex4/matlab/test2_11025.wav new file mode 100644 index 0000000000000000000000000000000000000000..04b882f23c78f86cdc6f6939b1d3ab256f70bbc5 Binary files /dev/null and b/Ex4/matlab/test2_11025.wav differ diff --git a/Ex4/matlab/test_signal.wav b/Ex4/matlab/test_signal.wav new file mode 100644 index 0000000000000000000000000000000000000000..51a5040ccd128813e34dbea06c96952f494151c2 Binary files /dev/null and b/Ex4/matlab/test_signal.wav differ diff --git a/Ex4/rundot.bat b/Ex4/rundot.bat new file mode 100644 index 0000000000000000000000000000000000000000..30e8d91718c1968819c6aab7a662145c94bc1e14 --- /dev/null +++ b/Ex4/rundot.bat @@ -0,0 +1,4 @@ +path = "D:\Program Files (x86)\Graphviz\bin";%path% + +del Ex4_task2.gif +dot -Tgif Ex4_task2.dot -oEx4_task2.gif