[my504 2]

18 minutes ago, dsl2640free said:

кусок кода, детектирующий гармоники в составе сигнала

Насколько я понял - это обычный ДПФ...

29 minutes ago, ViKo said:

Сюда ходи.

Благодарю. Вроде бы неплохо изложено в интересующем меня месте. Попробую еще раз зайти самостоятельно.

[my504 2]

Ну всё, разобрался. Всем спасибо.

Для dsPIC33 одна комплексная "бабочка" требует 20 машинных циклов из которых собственно mac-операции съедают 12, остальное уходит на адресацию исходных операндов и выгрузку результата. Оценочно, для dsPIC33E на 1024-точечный БПФ потребуется примерно 1,6 мс при максимальной скорости его работы (70 MIPS).

[my504 2]

Особенности указателей в DSP модуле dsPIC33 увеличили длину кода "бабочки" до  23 инструкций.

Получилось так:

; W10=pA
; W11=pB
; W8=pW
; A=A+B
; B=(A-B)W
BtflCmplx:
	mov		#0x7FFF, W4				; W4=0.99999
	clr		A, [W10]+=2, W6				; W6=Ar, [W10]->Ai
	clr		B, [W8]+=2, W5, [W11]+=2, W7		; W7=Br, [W11]->Bi, W5=Wr [W8]->Wi
;--- perform A
	mac		W4*W6, A, [W10]-=2, W6			; ACCA=+(1*Ar), W6=Ai, [W10]->Ar
	mac		W4*W7, A, [W11]-=2, W7			; ACCA=+(1*Br), W7=Bi, [W11]->Br
	sac.r		A, #-1, [W13++]				; Ar(new)=ACCA.rnd->buf
	mac		W4*W6, B, [W10]+=2, W6			; ACCB=+(1*Ai), W6=Ar, [W10]->Ai
	mac		W4*W7, B, [W11]+=2, W7			; ACCB=+(1*Bi), W7=Br, [W11]->Bi
	sac.r		B, #-1, [W13++]				; Ai(new)=ACCB.rnd->buf
;--- perform B
	mpy		W5*W6, A, [W10]-=2, W6			; ACCA=Wr*Ar, W6=Ai, [W10]->Ar
	msc		W5*W7, A, [W8]-=2, W5, [W11]-=2, W7	; ACCA=-(Wr*Br), W5=Wi, [W8]->Wr, W7=Bi, [W11]->Br
	msc		W5*W6, A				; ACCA=-(Wi*Ai)
	mac		W5*W7, A, [W8]+=2, W5			; ACCA=+(Wi*Bi), W5=Wr, [W8]->Wi
	sac.r		A, #-1, [W13++]				; Br(new)=ACCA.rnd->buf
	mpy		W5*W6, B, [W10]+=2, W6			; ACCB=Wr*Ai, W6=Ar, [W10]->Ai
	msc		W5*W7, B, [W8], W5, [W11]+=2, W7	; ACCB=-(Wr*Bi), W5=Wi, W7=Br, [W11]->Bi
	mac		W5*W6, B				; ACCB=+(Wi*Ar)
	msc		W5*W7, B				; ACCB=-(Wi*Br)
	sac.r		B, #-1, [W13]				; Bi(new)=ACCB.rnd->buf
	mov		[W13--], [W11--]			; save Bi
	mov		[W13--], [W11]				; save Br
	mov		[W13--], [W10--]			; save Ai
	mov		[W13], [W10]				; save Ar
	return

Итого - "бабочка" выполняется за 330 нс.

[dsl2640free 0]

Спойлер

 

Ниже приведен пример более сложной процедуры, которая использует DSP особенности dsPIC33E.

Спойлер

#include <dsp.h>;
#define FFT_BLOCK_LENGTH 1024 //Number of frequency points in the FFT
#define LOG2_BLOCK_LENGTH 10 //Number of "Butterfly" Stages in FFT processing should be related to FFT_BLOCK as in 2^9=512
#define AUDIO_FS 44100//Sampling frequency of audio signal captured by the mic
int16_t peakFrequencyBin;
uint16_t ix_MicADCbuff;
uint16_t peakFrequency;
extern const fractcomplex twiddleFactors[FFT_BLOCK_LENGTH/2]__attribute__ ((space(auto_psv), aligned (FFT_BLOCK_LENGTH*2)));
fractcomplex sigCmpx[FFT_BLOCK_LENGTH] __attribute__ ((section (".ydata, data, ymemory"), aligned (FFT_BLOCK_LENGTH * 2 *2)))={{0}};
void readMic(void)//Sample microphone input
{
ADC1_ChannelSelectSet(ADC1_AI_MIC);
ix_MicADCbuff=0;
for(ix_MicADCbuff=0;ix_MicADCbuff&lt;FFT_BLOCK_LENGTH;ix_MicADCbuff++)
{

//Insert delay subroutine here
ADC1_SamplingStop();
while(!ADC1_IsConversionComplete()){}
sigCmpx[ix_MicADCbuff].real = ADC1_Channel0ConversionResultGet();
sigCmpx[ix_MicADCbuff].imag = 0;
}
}

void signalFreq(void)//Detect the dominant frequency of the audio picked by the microphone
{
readMic();
FFTComplexIP (LOG2_BLOCK_LENGTH, &amp;sigCmpx[0], (fractcomplex *) __builtin_psvoffset(&amp;twiddleFactors[0]), (int) __builtin_psvpage(&amp;twiddleFactors[0]));// Perform FFT operation
BitReverseComplex (LOG2_BLOCK_LENGTH, &amp;sigCmpx[0]);// Store output samples in bit-reversed order of their addresses
SquareMagnitudeCplx(FFT_BLOCK_LENGTH, &amp;sigCmpx[0], &amp;sigCmpx[0].real);//Compute the square magnitude of the complex FFT output array so we have a Real output vector
VectorMax(FFT_BLOCK_LENGTH/2, &amp;sigCmpx[0].real, &amp;peakFrequencyBin);//Find the frequency Bin ( = index into the SigCmpx[] array) that has the largest energy
peakFrequency = peakFrequencyBin*(AUDIO_FS/FFT_BLOCK_LENGTH); //Compute the frequency (in Hz) of the largest spectral component
}

 

[my504 2]

Этот код использует не фичи dsPIC, а БИБЛИОТЕКУ, написанную кем то на АСМе. Потому что написать на голом Си работу с DSP ядром этой платформы совершенно невозможно. Полагаю, что и на аналогичных архитектурах DSP других производителей - тоже. Даже не представляю себе, как заставить компилятор управлять указателями и загрузкой по ним в инструкциях DSP...

[ilyamp 0]

Строго говоря, БПФ тут реализовано двумя строками (даже одной):

On 1/27/2020 at 3:50 PM, dsl2640free said:

FFTComplexIP (LOG2_BLOCK_LENGTH, &amp;sigCmpx[0], (fractcomplex *) __builtin_psvoffset(&amp;twiddleFactors[0]), (int) __builtin_psvpage(&amp;twiddleFactors[0]));// Perform FFT operation
BitReverseComplex (LOG2_BLOCK_LENGTH, &amp;sigCmpx[0]);// Store output samples in bit-reversed order of their addresses

сказать, как и насколько функция  FFTComplexIP (...) учитывает особенности dsPIC33, на мой взгляд, нельзя. 

если на асме прилично оптимизировать бабочку(правильные предвыборки и логичные смещения указателей), то дальше ускорить БПФ практически некуда

[DASM 0]

On 1/27/2020 at 5:58 PM, my504 said:

Этот код использует не фичи dsPIC, а БИБЛИОТЕКУ, написанную кем то на АСМе. Потому что написать на голом Си работу с DSP ядром этой платформы совершенно невозможно. Полагаю, что и на аналогичных архитектурах DSP других производителей - тоже. Даже не представляю себе, как заставить компилятор управлять указателями и загрузкой по ним в инструкциях DSP...

 Выделяет характерные паттерны и в путь. Когда писал (немного) для Shark - тоже удивлялся, он и инструкции правильно распараллеливал, и аппаратные циклы использовал. Тоже самое сейчас вижу код под Neon создаётся- править после практически нечего. А скоро думаю AI направят на все это дело и программистам останется только ТЗ грамотно писать.