Tuesday, August 23, 2022

DFT_grid.m

Pitch detection by DFT grid
 
Acknowledgements

1. In the late 80's and early 90's, I studied Byzantine music in the byz music school of Piraeus (Greece). There, I learned the basics of byz music solfège and I decided to create a project to help me learn to sing.
2. In 90's and in 2004 I was studing DSP books in Eugenidis Foundation Library. There, I learned DSP basics.
3. In 2016 I studied Mr. Richard Lyons, a mentor for DSP specialists, book: Understanding Digital Signal Processing, and learned DFT.
4. In 2017 I asked at dsp.stackexchange question: Can I test any frequency with DFT? Marcus Müller, also mentor for DSP specialists, answered *yes*. So I learned can perform DFT grid.
5. Mr. Jeffrey Clark, at matlab suggested me 16384 for second argument in function fft.
 
Introduction

1. I'm homeless and can't implement this code in DSP kit. I offer it for free. It detects fundamental frequency (fundFreq), or pitch, or 1st harmonic, using DFT grid about fundFreq detected by FFT or any other method detecting pitch. My method is the condition user give pitch between one tone (9/8) lower than note's frequency and one tone higher than it, even a student of music can do. Then program, by matlab's abs(fft()), find indices up to 1.2 note's frequency. 1st harmonic's abs is the maximum.
2. Indeed it should be embedded in a project executing, in part, in DSP kit.
3. Project should have two parts. One, executing in DSP kit
performs DFT grid, and other, executing in PC's processor, performs the graphics.
4. Graphics part represents musical note as rectangle with height representing note's interval from base by some unit e.g. cent of semitone. Wide represents note's duration (temp). Color represents sound's amplitude, volume. Rectangle moves to left with speed related with wide. At the left side of window a pointer moves up and down. Its point's height from base, represents voice's pitch. For get an idea see Metrw-pitch graphics. 
5. You can create vertical scale with zero at its middle, upper half representing one semitone, and lower also one semitone, just a micrometer. Zero represents the note rectangle's upper edge. A speed regulator is a must, as well a horizontal scroll bar for user can go back.
6. To measure the volume of the sound, we need to add a snippet. Also we need a third scale for sound's volume.
 
Accuracy error
 
1. Singers, chanters, violists have to give even the smallest musical interval, diesis. But fundFreq, detected by FFT, has accuracy error about a diesis, so projects detecting fundFreq by FFT do not satisfy musicians. Similar errors occur in other methods detecting pitch. Present code detects funbdFreq with accuracy better than 0.1 diesis. So when this code is embedded in a project executing in DSP kit, user can see in real time 8 times per second his voice's pitch error and regulate voice's pitch for give diesis. Question is whether MATLAB's code executing in dsp kit, retains same accuracy as in _MATLAB online R2020b. See answer from Mr. Ben to my question at dsp.stackexchange.
2. Error depends on the ratio Segm/SampFreq. As greater the ratio as higher the accuracy. But as lower the pitch's iteration rate in graphics.
3. Error depends also on fundFreq. As lower fundFreq as greater accuracy error. Reason for both is obvious: As more FundFreq's periods in Segm, as better accuracy. Moreover it depends, to a lesser degree, on amplitude, phase, and upper harmonics.
 
General instructions

1. For avoid analog filter which distorts signal, I suggest sampling frequency 96000 samples/sec. Some noise reduction snippet is necessary. After that a
2. Low pass FIR filter with transition band 7000-8000 khz.
3. Choosing 1 sample every 6, SampFreq is reduced to16000 samples/sec.
4. Programmer can change numbers.
5. For a conical presentation of FFT amplitudes see answers to my question at dsp.stackexchange.
 
 
The code    

format default;
SampFreq = 16000;
Segm = 1:2048;
baseFreq = 256;
noteFreq = baseFreq*100/81; % up to 2666, for 3rd harmonic < 8000
x = 1.1; % 8/9 < x < 9/8
voicePitch = x * noteFreq;
% Generating signal
FirstHarmAngles = voicePitch*2*pi/SampFreq*Segm+1.9*pi;
SinFirstHarmAngles = sin(FirstHarmAngles);
SecondHarmAngles = voicePitch*2*2*pi/SampFreq*Segm+2.9*pi;
SinSecondHarmAngles = sin(SecondHarmAngles);
ThirdHarmAngles = voicePitch*3*2*pi/SampFreq*Segm+0.3*pi;
SinThirdHarmAngles = sin(ThirdHarmAngles);
Xn=170000*SinFirstHarmAngles+220000*...
SinSecondHarmAngles+150000*...
SinThirdHarmAngles;
% Finding FFT pitch
magn = abs(fft(Xn, 16384));
OnePointTwo = round(noteFreq * 1.2);
[maxMagn, FFTpitch] = max(magn(1:OnePointTwo));
% Performing DFT grid
LowSpanLimit = 0.8*FFTpitch;
HighSpanLimit = 1.2*FFTpitch;
PartsOfSpan = 1001;
TestFreqs = logspace(log10(LowSpanLimit), log10(HighSpanLimit), PartsOfSpan);
MagnSqrd = ones(1, PartsOfSpan);
for f = 1:PartsOfSpan
Angles = TestFreqs(f)*2*pi/SampFreq*Segm;
XnCos = sum(Xn.*cos(Angles));
XnSin = sum(Xn.*-sin(Angles));
MagnSqrd(f) = XnCos.^2+XnSin.^2;
end
% Finding GRID pitch
[maxMagnSqrd, maxMagnSqrdIndex] = max(MagnSqrd);
GRIDpitch = TestFreqs(maxMagnSqrdIndex);
voiceRatioInterval = GRIDpitch/baseFreq;
% Following should be incorporated in graphics part.
mainUnit = 9/8; % for cents of semitone it is 2
tempUnit = 100; % for cents of semitone it is 1200
noteRatioInterval = noteFreq / baseFreq;
tempNOTEpitch = log(noteRatioInterval)/log(mainUnit)*tempUnit;
tempGRIDpitch = log(voiceRatioInterval)/log(mainUnit)*tempUnit;
% Display results
noteFreq
voicePitch
FFTpitch
GRIDpitch
tempGRIDpitch
tempNOTEpitch