PyTTaMaster / PyTTa

Python in Technical Acoustics and Vibration
MIT License
92 stars 32 forks source link

Using loopback in FRF measurement #31

Closed labodezao closed 3 years ago

labodezao commented 3 years ago

Hi,

in the reaserch gate paper : https://www.researchgate.net/publication/336207488_PyTTa_Open_source_toolbox_for_acoustic_measurement_and_signal_processing

It is said that it's possible to use a loopback, certainly to avoid shift of frequency. But I did not find any test / exemple in the docs / forums to explain how to implementate this in the code.

I am actually making a FRF from vibrations, input is a electromagnéctic exiter (which works at Frequency /2 due to magnetic pole ) and I have a accelerometer to measure the response. I have a shift of frequency about 10Hz right now (I compare with a signal generator and the Pyyta results)

I wonder also if it's possible to make a FRF from impulse excitation.

Thanks a lot

My code is inspired by frf exemple below :

import pytta

if __name__ == '__main__':

    device = pytta.get_device_from_user()

    sweep = pytta.generate.sweep(fftDegree=20,samplingRate=15000, freqMin=39.4, freqMax=42)  # Generate sine sweep SignalObj
    #sweep = pytta.generate.sweep(fftDegree=20, samplingRate=15000, freqMin=41, freqMax=43,startMargin=0.05, stopMargin=0.05,method='linear')  # Generate sine sweep SignalObj

    measurementParams = {
        'excitation': sweep,  # Passes the sweep signal as excitation signal
        'samplingRate': 15000,  # Frequency of sampling
        'freqMin': 20,  # Minimum frequency of interest NOT WORKING
        'freqMax': 20000,  # Maximum frequency of interest NOT WORKING
        'device': [1,5],  # Device number provided at runtime
        'inChannels': [1,2],  # List of hardware channels to be used
        'outChannels': [1],  # List of hardware channels to be used
        'comment': 'Testing; 1, 2.'  # just a comentary
    }

    ms = pytta.generate.measurement('frf',  # Generates the configuration for an impulse response measurement
                                    **measurementParams)

    #%% Run the measurement and process data, return an ImpulseResponse object
    med1 = ms.run()

    #%% Shows the data
    med1.plot_time()
    med1.plot_freq(smooth=True)

    #%% Save the measurement
    path = 'data\\'
    name = 'myir'
    filename = f'{path}{name}'

    # Save PyTTaObj as HDF5 file
    pytta.save(filename, med1)

    # Export wave file of SignalObj
    pytta.write_wav(filename+'.wav', med1.IR)
Chum4k3r commented 3 years ago

Hi @labodezao, thank you for issuing this here.

The loopback mentioned in the paper is actually a FRF of the ADC/DAC hardware used on the measurements.

To acquire it you must connect the output of the DAC to the input of the ADC and run a measurement just like you did in your code.


# Adapted from the code you provided.

import pytta

if __name__ == '__main__':

    device = pytta.get_device_from_user()

    sweep = pytta.generate.sweep(fftDegree=20,samplingRate=15000, freqMin=39.4, freqMax=42)  # Generate sine sweep SignalObj
    #sweep = pytta.generate.sweep(fftDegree=20, samplingRate=15000, freqMin=41, freqMax=43,startMargin=0.05, stopMargin=0.05,method='linear')  # Generate sine sweep SignalObj

    measurementParams = {
        'excitation': sweep,  # Passes the sweep signal as excitation signal
        'samplingRate': 15000,  # Frequency of sampling
        'freqMin': 20,  # Minimum frequency of interest NOT WORKING
        'freqMax': 20000,  # Maximum frequency of interest NOT WORKING
        'device': device,  # Device number provided at runtime
        'inChannels': [1,2],  # List of hardware channels to be used
        'outChannels': [1],  # List of hardware channels to be used
        'comment': 'Testing; 1, 2.'  # just a comentary
    }

    ms = pytta.generate.measurement('frf',  # Generates the configuration for an impulse response measurement
                                    **measurementParams)

    input('At this point, make the direct DAC-to-ADC connection. Press Enter to continue')
    loopbackSignal = ms.run()

    #%% Run the measurement and process data, return an ImpulseResponse object
    input('Now connect your measurement setup as you would normally. Press Enter to continue')
    med1 = ms.run()

    # To apply the loopback correction, first you must divide the med1 by the loopbackSignal
    corrMed = med1.IR / loopbackSignal.IR  # recorded signal with loopback correction

    #%% Shows the data
    corrMed.plot_time()
    corrMed.plot_freq(smooth=True)

    #%% Save the measurement
    path = 'data\\'
    name = 'myir'
    filename = f'{path}{name}'

    # Save PyTTaObj as HDF5 file
    pytta.save(filename, med1, loopbackSignal, corrMed)

    # Export wave file of SignalObj
    pytta.write_wav(filename+'.wav', corrMed.IR)
labodezao commented 3 years ago

I have another question about my code. As you can see I use a sampling rate of 1500 Hz because my system is mechanical and I need the measurement to be slower that the interial constant time of my system. Do I need to have the exactly same sampling rate for excitation and measurement signals ? Do slower sampling rate could generate shift frequency error in the frf curve ? For structural vibration analysis, what is the good / best sampling frequency ?

Thanks a lot !

Chum4k3r commented 3 years ago

Whenever you are operating signals, they must be on the same sampling rate.

The proper sampling rate is the one your soundcard accepts. Right now PyTTa does not provide any means of finding out which sampling rates the soundcard can work with. There are some with fixed sampling rate, say 44100 or 48000 Hz, and there are others that can deal with ranges of sampling rate values, as long as it stays fixed during runtime, say Fs <= 48000 Hz, so you can use unusual values such as 1500 Hz.

I suggest you try some common sampling rate as a test case, like 44100 Hz. Then try 8000 Hz, if it works, go for 2000 Hz.

labodezao commented 3 years ago

I use a Behinger uCONTROL https://www.thomann.de/fr/behringer_ucontrol_uca_222.htm in my windows 10 driver I put 48.000 KHz and same in pytta

Regarding loopback, i have a card with 2 inputs 2 outputs. Is it possible to make loopback and measurement at the same time ou is it better to male the loopback over the input / output used in measurements ?

Thanks a lot !

labodezao commented 3 years ago

This is my result... Weird because the resonnance is always 50 Hz ... Like the electrical french oscilatting frequency ...

labodezao commented 3 years ago

It seems that the signal is quite close to Zero .... image

labodezao commented 3 years ago

Seems that the problem is that Pytta doesn't select the proper input number in my souncard...

Do you have a test script to record the input to see if everything in right in the soundevices ?

I have several option ....

2 Microphone (2- USB Audio CODEC , MME (2 in, 0 out)

7 Haut-parleurs (2- USB Audio COD, MME (0 in, 2 out)

10 Microphone (2- USB Audio CODEC ), Windows DirectSound (2 in, 0 out)

15 Haut-parleurs (2- USB Audio CODEC ), Windows DirectSound (0 in, 2 out)

18 Haut-parleurs (2- USB Audio CODEC ), Windows WASAPI (0 in, 2 out)

21 Microphone (2- USB Audio CODEC ), Windows WASAPI (1 in, 0 out)

28 Haut-parleurs (USB Audio CODEC), Windows WDM-KS (0 in, 2 out)

29 Microphone (USB Audio CODEC), Windows WDM-KS (2 in, 0 out)

Which driver is best ?

My signal is on the Left input of my soundcard and I get high level signal in audacity :)

`import pytta

if name == 'main':

device = pytta.get_device_from_user()

sweep = pytta.generate.sweep(fftDegree=20,samplingRate=44100, freqMin=39, freqMax=42)  # Generate sine sweep SignalObj
#sweep = pytta.generate.sweep(fftDegree=20, samplingRate=15000, freqMin=41, freqMax=43,startMargin=0.05, stopMargin=0.05,method='linear')  # Generate sine sweep SignalObj

measurementParams = {
    'excitation': sweep,  # Passes the sweep signal as excitation signal
    'samplingRate': 44100,  # Frequency of sampling
    'freqMin': 20,  # Minimum frequency of interest NOT WORKING
    'freqMax': 50,  # Maximum frequency of interest NOT WORKING
    'device': [2,7],  # Device number provided at runtime
    'inChannels': [1],  # List of hardware channels to be used
    'outChannels': [1],  # List of hardware channels to be used
    'comment': 'Testing; 1, 2.'  # just a comentary
}

ms = pytta.generate.measurement('frf',  # Generates the configuration for an impulse response measurement
                                **measurementParams)

#%% Run the measurement and process data, return an ImpulseResponse object
med1 = ms.run()

#%% Shows the data
med1.plot_time()
med1.plot_freq()

#%% Save the measurement
path = 'data\\'
name = 'myir'
filename = f'{path}{name}'

# Save PyTTaObj as HDF5 file
pytta.save(filename, med1)

# Export wave file of SignalObj
pytta.write_wav(filename+'.wav', med1.IR)

`

labodezao commented 3 years ago

I did not find any solution yet (it seems that the value recorded of pytta is close to 0) I join here the differents value grabbed from my osciloscope and Pytta if you see something I did not see ...

Signal CH1 : burst

Signal CH (blue) Accelerometer

The Pytta signal is also joined to this message... and this is not looking really the same that what I get on my osciloscope :)

osc 2 sigs oscilos with values temp pytta

Chum4k3r commented 3 years ago

Let's go bottom up:

The Pytta signal is also joined to this message... and this is not looking really the same that what I get on my osciloscope :)

There is a difference in what you read in a Osciloscope and what a FRF measurement delivers. The Osciloscope shows the record signal, while the FRF shows you the signal resultant of dividing the record by the excitation. This is also the reason why it has such low amplitude.

Seems that the problem is that Pytta doesn't select the proper input number in my souncard...

Do you have a test script to record the input to see if everything in right in the soundevices ?

I have just the same soundcard. And it is hell to use it on windows hahahaha the lack of a proper driver make it a bit tricky to work with on stuff like these, but it is not impossible. To find the best device you'll need to go one on one and try all driver IDs that represent the 'USB Audio Codec'. Usually the DirectSound works better for me, but it varies. Try easing the load on your computer, not many software running at the same time you are running your measurements.

This is my result... Weird because the resonnance is always 50 Hz ... Like the electrical french oscilatting frequency ...

If you are reading a peak at the same frequency of your electrical energy source of, say, your house or laboratory, it is a common issue on audio and electronics. It is not a frequency shift, it is a signal leakage. In this case, the electric signal of some power cord is leaking through your equipment electronics. Try recording without energy cables (of course only on equipment which has batteries). Like if you are using a laptop, try removing the power cord while you run your measurement.

Regarding loopback, i have a card with 2 inputs 2 outputs. Is it possible to make loopback and measurement at the same time ou is it better to male the loopback over the input / output used in measurements ?

It only makes sense to do the loopback on the same channels you are using to run your measurements.

labodezao commented 3 years ago

Hi ! I have (few ?) more question :) 👍

1) Is it possible to plot the oscilloscope like signals in Pytta ? It could be really nice to see if the measurement looks good :)

2) On my oscilloscope, I see that my response signal have a hudge continus component ( mean value = 770 mV) when the alternative component is only 128 mV ... So I don't know if the continus component of the signal is filtered by the sound card but can this could make troubles ?

3) Compared to the excitation the response signal is so verry low, is it possibler to put an attenuation before making the division on the excitation in order to have similar scale for the signals magnitude ? Because if I divide a low 128mV signal by my huge

4) Is it possible to perform a FRF / modal analysis with a step / impulse excitation with Pytta ?

On the software Room eq Wizard, I did the same test and there is my results : image

I expect my results looking like this : https://github.com/openmodal/pyFRF/blob/master/Showcase%20pyFRF.ipynb

Best and Thanks

Chum4k3r commented 3 years ago

As I've said, the osciloscope shows the recorded signal, so you can retrieve something similar from med1 object if you access the med1.record attribute, which is a SignalObj of the recorded signal. You can visualize it just like the IR.

3. Compared to the excitation the response signal is so verry low, is it possibler to put an attenuation before making the division on the excitation in order to have similar scale for the signals magnitude ? Because if I divide a low 128mV signal by my huge

Well, mathematically you can rescale it by any factor. All you have to do is multiply either the recorded signal or the resulting IR. med1.record.timeSignal = med1.record.timeSignal * 4 for example. What it will lose is some sort of physical representation of how all your measurement system (soundcard, exciter, accelerometer, etc) alters the magnitude of the signal you are playing.

4. Is it possible to perform a FRF / modal analysis with a step / impulse excitation with Pytta ?

The exponential sine sweep waveform you are generating is an extremely precise solution for this modal analysis problem. Right now we don't have a sine step generator, although it is a great idea. The impulse you can create with pytta.generate.impulse() but it is more likely to be a theoretical signal rather than a useful signal. Feel free to try it out, though!

What I do believe you are looking for is actually the .plot_freq() method, like in med1.IR.plot_freq() so you can see the Frequency Response, and not the Impulsive Response.

labodezao commented 3 years ago

As I've said, the osciloscope shows the recorded signal, so you can retrieve something similar from med1 object if you access the med1.record attribute, which is a SignalObj of the recorded signal. You can visualize it just like the IR.

Nice You :) Good to know, That was not that obvious in documentation !

med1.record.timeSignal = med1.record.timeSignal * 4 for example. What it will lose is some sort of physical representation of how all your measurement system (soundcard, exciter, accelerometer, etc) alters the magnitude of the signal you are playing.

Hmm but If I use a preamp for amplify the accelerometer signal before the soundcard as it is made habitually, it will alters also the physical representation of the signal no ? I will try to normalize my results So :)

Is it possible to make the sweep from higher frequency to lower ? like start a 1000 Hz and finish at 100 Hz ?

Right now we don't have a sine step generator, although it is a great idea.

I think that I need to make 3 kinds of measuements with my setup and I know from experience (I am reasearcher and also accordion maker) that the self oscillations threshold have a kind of hysteresis. I am sure that a sine step generator could be great for transients. What could be the application of a step sine generator ?

What I do believe you are looking for is actually the .plot_freq() method, like in med1.IR.plot_freq() so you can see the Frequency Response, and not the Impulsive Response.

What is the difference between IR and FRF ? I mean the Impulsive response ans the Frequency response ? I will be happy to write a paper based on those experiements, but actually I am a 'freelancer' scientist so if you are intersted by a collaboration and If that could help you I will be pleased to :)

Thanks again !

labodezao commented 3 years ago

Hi ! I can't find out where the recorded signals object are stored for FRF measurement !

med1.record.timeSignal = med1.record.timeSignal * 4 for example. What it will lose is some sort of physical representation of how all your measurement system (soundcard, exciter, accelerometer, etc) alters the magnitude of the signal you are playing.

When I type med1.record.timeSignal I have an error that

If I plot : med1.irSignal.plot_time() med1.systemSignal.plot_time() med1.tfSignal.plot_time()

I get exaclty the same curves... So I don't know where are stored the excitation and recorded signal and I don't know what the différences between rSignal., systemSignal. and .tfSignal.

labodezao commented 3 years ago

Hi, I made a return in the run() FRF fonction and I was able to have the record SignalObject function. Maybe it could be nice to have the possibility to access it without modify the code. I don'yt know Github so I don't know if what I made is good, but it works... Now I think the results are really weird : image

Here it is my sweep function and my record function (That was clipping but I could not see the range of it before) And there is my results of the IR computing made by PyTTa : image

And the time signal looks really weird and the frequency response have no spike wheras I have a system with big resonances and a Q factor at least 25 ! The resonnance bandwith is really narrow ( maybe 0.1 Hz) Maybe I should put a nfft geater than I have (= 22)

Waiting from your answear, I wish you the Best

labodezao commented 3 years ago

Hi, I wh=ish you a good monday. Antonin Novak, A French reasercher from LAUM from where I am graduate wrote this paper : https://ant-novak.com/pages/sss/ This is about the improvement in the way to generate sweep sine exponential function and to perform the deconvolution in frequency domain or in time.

I wonder if Pytta use this kind of deconvolution method, but I think I will try to see what is the results with this methods and will let you know if I get a Pike in my FRF.

I will be happy to have your feedback on my last question. Have a good Monday !

Chum4k3r commented 3 years ago

We do have Angelo Farina correction implemented on deconvolution.

But we have not made any code for non linear systems. Our first study object was rooms transfer functions, which behave as time invariant linear systems.

If you have any idea or intentions on bringing on this non linear system evaluation method, I do believe will be a great contribution, just as the other you suggested.

Feel free to make a PR.

mtslazarin commented 3 years ago

I get exaclty the same curves... So I don't know where are stored the excitation and recorded signal and I don't know what the différences between rSignal., systemSignal. and .tfSignal.

Please, read the documentation. https://pytta.readthedocs.io/pyttadoc/classes/signal.html

I made a return in the run() FRF fonction and I was able to have the record SignalObject function.

You can always use a PlayRecMeasure instead a FRF measurement and calculate the impulsive response by yourself through the ImpulsiveResponse class.

And the time signal looks really weird and the frequency response have no spike wheras I have a system with big resonances and a Q factor at least 25 !

The time signal is weird because you have too much noise. You should set the minFreq and maxFreq parameters in the impulsive response calculation according to the analysis you want to do, then the Kirkeby regularization will take care of the noise. You should review the signal processing concepts involving what you are trying to do.