Jump to content

    

Спектральный анализ Фурье

18 minutes ago, dsl2640free said:

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

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

29 minutes ago, ViKo said:

Сюда ходи.

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

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

Особенности указателей в 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 нс.

Share this post


Link to post
Share on other sites
Спойлер

 

Ниже приведен пример более сложной процедуры, которая использует 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
}

 

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

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

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, на мой взгляд, нельзя. 

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

Share this post


Link to post
Share on other sites
On 1/27/2020 at 5:58 PM, my504 said:

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

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

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now