Fourier Transform of sampled reflectance measurements Goal:
Calculate the complex index of refraction ($\hat n = n +jk)$ from reflectance measurements.

Data:

Reflectance measurements for various materials were taken with a FTIR (Nicolet™ iS50 FTIR, from Thermofisher). The spectral range of the reflected measurements are from 400nm - 15um, they are reported in both wavenumber and wavelength (relationship: $\hat\nu = \frac{1}{\lambda} [cm^{-1}]) $. Note, that  the full wavelength range was captured with three different detectors.

To simplify the matter (and give access to data) I have been using reflectance measurements form the following NASA database https://speclib.jpl.nasa.gov/library. I have chosen SiO2 from the database (2um-15um) which can be navigated to by selecting:

Select Spectral Type-> Minerals

Filter by class-> Silicate

Wavelength-> VSWIR+TIR

Quartz SiO_2

Brief Theory:

Here's a short summary: https://shimadzu.com.au/sites/default/files/Appl_FTIR_Polymer_specular_reflectance_055_en.pdf

The imaginary part of an analytical signal can be found from the real part alone through the Hilbert transform, e.g. $\tilde {x_c}(t) = x_r(t) +j x_i(t)$ meaning $x_i(t) = HT{x_r(t)} $. https://epdf.pub/hilbert-transforms-volume-1-encyclopedia-of-mathematics-and-its-applications.html

Changing mindset to the optical domain, and the question at hand, noting that most optical literature calls the Hilbert transform the Kramer Kronig. The FTIR measures spectral reflectance (i.e. Intensity).

$$\tilde r = re^{j\phi} = \sqrt{R}e^{j\phi}$$

$$R = |\tilde r|^2$$

Where r is the reflectivity
R is the reflectance (This is what is measurable, intensity)
$\phi$ is the phase change at the surface caused by the absorption of the material.

Through the Fresnel equations the complex refractive index can be calculated from the below equations.

Dispersion:
$$n(\nu) = \frac{1-R(\nu)}{1 + R(\nu) - 2\sqrt{R(\nu)}cos(\phi(\nu))}$$

Absorptive index:
$$k(\nu) = \frac{-2\sqrt{R(\nu)}sin(\phi(\nu))}{1 + R(\nu) - 2\sqrt{R(\nu)}cos(\phi(\nu))}$$

So, the goal is to take the measured reflectance data and calculate the phase through the Hilbert Transform, as shown below.

$$\phi(\nu_g) = \frac{2\nu_g}{\pi} \int_0^\infty \frac{ln\sqrt{R(\nu)}}{\nu^2 - \nu_g^2}$$

where $\nu$ is the wavenumber

Since I have discrete data (over a finite range) the Hilbert transform can't be calculated directly. It's usually calculated through Maclaurin's method or the double Fourier Transform (Discrete Fourier Transform (DFT) in this case).

The approximation of the Hilbert transform using the double FT is given shown in the below equation.

$$\phi(\nu_g) = 4 \int_0^\infty cos(2\pi \nu_gt)dt \int_0^\infty ln\sqrt{R(\nu)}sin(2\pi \nu_gt)dv$$

My Question:

My question pertains to the Fourier transform (FT) of the sampled data. All the equations I see on line show the FT of evenly spaced sample data with respect to "time", which I don't know. My data is evenly sampled with respect to wavenumber $\nu$, can I still use the FT but replace the sampling time (1/$\Delta$t) with (1/$\Delta\nu$)? If yes, is there anything extra I need to consider with the respect to the reluctance data.

The last equation showing the double Fourier transform has a intergal with respect to time, but I don't have any sampling time information, is there a way around this, or some conversion?

I found the following relationships online but can't quite see how to apply them.

$$f = \frac{1}{t} ....||||.... \nu = \frac{1}{\lambda}$$

$$w = 2\pi f ....||||.... k = \frac{2\pi}{\lambda}$$

This elaborates on my comment and my linked FTIR example. Once the FTIR spectrum is obtained, the phase spectrum can be obtained via the Hilbert Transform of the spectrum and that can be performed as a convolution. The Hilbert Transform can be performed by various software packages. In Igor Pro (v. 6.3), the Help text on Hilbert Transforms (page 723 of the Manual), is simply:

In the linked example, the IR spectrum of ethyl acrylate was obtained via a Mueller optical calculus simulation. For present illustrative purposes, the spectrum was adjusted in two ways. First, everything below 800 wavenumbers was set to zero. This is typically where the noise is worst because IR detectors are relatively insensitive and noisy compared with detectors in higher energy ranges. Second, the spectrum was zero filled from 4001 to 4095 wavenumbers, in order to achieve a total of 4096 wavenumber values. This spectrum is shown below:

Performing the Hilbert Transform, as per the Igor Pro command in the first figure, then yields the phase spectrum:

Maybe this helps and maybe it doesn't. But it points to a way forward: find a reflectance spectrum in the refereed literature and the dispersion and absorption spectra that were computed from it. This has to be something about which there is no doubt, i.e., it will be the standard. Then try using software to compute the Hilbert transform and see if the result agrees with the standard. Lastly, if you still want to do so, then pursue the pathway you have outlined in your question text. Hope this helps a bit!


How can you 'perfectly' reconstruct phase for transients in a modified audio absolute spectrogram? I am currently working on the absolute spectrogram of an input audio signal. I am using non-negative matrix factorization to separate out the rhythmic components (bass drum, snare drum and hi-hat) of the file and shift around the activations of these components to get a modified spectrogram.

Now, I want to reconstruct the phase for this modified spectrogram so that I can take an inverse FFT and recreate the audio signal in time domain. I am currently using the Griffin Lim and the Phase Gradient Heap Integration methods from this phase retrieval toolbox. Unfortunately, this doesn't produce good results. Maybe it's with the arguments I'm using. I'm currently using a window size of 2048 with a hop of 512.

1, 2 - These questions were previously posted where the idea of reconstruction of audio from absolute spectrogram was tackled. This question pertains to how well one can reconstruct the audio signal back in the time domain and what would be the best way to do so when the harmonic content would be the same but the transients would be time-shifted. Also, how would you objectively or subjectively evaluate such a system formally?

Thanks!

Real Time Iterative Spectrum Inversion would be a way to tackle that problem. Have a look at this paper:

IMPROVING RTISI PHASE ESTIMATION WITH ENERGY ORDER AND PHASE
UNWRAPPING


How to estimate and compensate for doppler shift in wireless signals? I was wondering what good method(s) exist for the estimation (and subsequent compensation) of doppler shift for transmitted signals, be they acoustic or RF, in the context of comms.

The question: Specifically, if the degree of doppler shift is varying over the duration of a packet, how best to estimate (track it?) and then compensate for it. Assume we have a trainer sequence. You may also assume the signal passband BW is on the order of its carrier. (e.g, if the passband signal exists from 2500-7500 Hz, its BW is 5000Hz, as is its carrier.)

Some additional background for context:

One method I found during my research:

Since I have a trainer sequence and know its frequency, first I estimate its received      frequency.
Next, I resample the entire packet by a ratio that is related to the speed of the wave in the medium, my known transmitted frequency, and my new estimated doppler shifted frequency.
This works well in simulations, but the weak points are that the frequency estimation must be very accurate, and it also assumes the doppler shift is not changing during the duration of the packet.

Are there any other methods that can be utilized to solve the problem when the doppler is changing during the duration of the packet? What is the opinion of the above aforementioned method?

Many thanks!

This is a very common communications problem. Look in a textbook for "frequency synchronization"; whole books have been written on these and related topics. The technique that you would choose is a function of the specifics of your system. There are two common sources of frequency offset:

Differences in frequency between the reference oscillator at the transmitter and receiver. This error is usually small, depending upon the precision of the timebases available, and can be mitigated at some cost. Cheap crystal oscillators typically achieve 50 parts per million of error or better (although this does drift as the crystal ages). If you have a larger budget, you can use something like a Rubidium standard, which provides ~1 part per trillion of frequency error. A cheaper and increasingly common approach is to use a GPS receiver that has precision frequency output (typically 10 MHz). The highly accurate timebase available from the GPS constellation can be used to train the reference precisely on-frequency.
Physical dynamics effects between the transmitter and receiver. One notable example in which this comes into play are for satellite communications applications (especially in lower orbits), where the satellite is moving (and accelerating) very rapidly relative to any observer on the Earth. The high radial velocity of the satellite toward the receiver will cause a Doppler shift, and any change in that velocity caused by its orbit will cause that shift to change over time. In applications where you have these sort of dynamics, you typically can't mitigate it much, so you're left with building a receiver that tolerates the effects.

So how does a receiver synchronize with the transmitter in these cases?

One common approach useful for phase- or frequency-modulated signals is to use a phase-locked loop. Design of PLLs is a complex topic in itself, but in essence, they are feedback systems can be used to acquire and track phase and frequency offset while your receiver is in operation. If you only need frequency synchronization, then you can use a frequency-locked loop instead; while they will not provide phase synchronization for you, they often have better acquisition properties.
As an alternative to a feedback loop, there are also feedforward approaches to estimating frequency or phase offset. One feedforward approach would take advantage of your training sequence to estimate the frequency error based on how the phase offset changes over the course of the sequence. If the frequency offset changes as time goes on, however, you would need to repeat the estimation procedure to allow your receiver to catch up.
Another technique is to design your system to be robust to (reasonably small) frequency offsets. Differentially-encoded phase modulation is an example of this (although frequency error will show up as a phase offset after differential decoding, which must be handled). Frequency-modulated waveforms such as FSK also have some level of resistance to frequency offset, as long as the offset is small relative to the amount of frequency deviation used by the transmitter.

This very short summary really just scratches the surface of some of the more well-known approaches. Synchronization can be a difficult problem to solve in practical ways, and there has been much research over the years on different ways to do it. It's going to depend upon exactly how your system is structured, and one very important variable: the target SNR. There isn't a single "correct" answer. I will make one textbook recommendation; while it is very expensive, "Synchronization Techniques for Digital Receivers" by Mengali is a comprehensive text on timing, phase, and frequency synchronization.


How to allign audio signals after DTW? Hello I am not a domain expert in signal processing but I need to align two audio signals.

I have seen the following page
https://librosa.github.io/librosa_gallery/auto_examples/plot_music_sync

I have implemented DTW and so many other pages. Most of these are using DTW but the problem is that none of them are showing WHAT TO DO AFTER I PERFORM DTW!. DTW gives you a cost matrix and wrapping path (D, wp) but nowhere it is shown how to use these two to align the actual signals.

I have come to a github repo
https://github.com/chausies/voice_transmogrifier/blob/master/sync_audio.py

But most of the parts of the code seems difficult to understand. Like why he used variable_phase_vocoder or what close_points is doing. Why to use PchipInterpolator ? In fact why interpolation is done at all in audio?

Please help me in the project. I am new to audio processing.

Thanks

none of them are showing WHAT TO DO AFTER I PERFORM DTW!.

There are two quotes, at least in the page that you have linked that hint at what you do next, these are:

"We assume that you are familiar with the algorithm and focus on the application. Further information about the algorithm can be found in the literature..."

The linked reference is not a bad start for a book on audio signal processing. The second quote is towards the end:

"...you could apply time scale modification algorithms, e.g. speed up the slower signal to the tempo of the faster one."

And this provides a way towards the answer to your secondary question too.

The "cost matrix" that Dynamic Time Warping returns is discrete with units of samples (e.g. match sample-to-sample), or time.

With reference to this diagram from the Librosa documentation:

You can see that this "picture" has two axes: The $x$ and $y$ axis. On $x$ axis, are the time instances that samples occur on one waveform. On the $y$ axis, are the time instances that the same samples occur on the second waveform.

If the two waveforms were identical, their samples would be played at the same instances and the cost matrix would describe a straight line which means that sample 0 on $x$ is found on sample 0 on $y$, sample 1 on $x$ is found on 1 in $y$ and so on until sample $N$.

But, most of the times you apply DTW, you are interested in aligning the two waveforms.

In that case, the "cost matrix" can be interpreted as a speed-up or slow-down factor, so that the samples of waveform $x$ catch-up with those of waveform $y$.

  ...why interpolation is done at all in audio?

When you think of the cost matrix as modifying the speed by which samples are reproduced, you start thinking about how to speed up or slow down a waveform.

Especially in audio, there are many different methods by which you can achieve this. But the easiest way to think about it is as if changing the speed by which your "sample pointer" rolls through the waveform.

The most trivial example is this: At the waveform's $Fs$ (Sampling frequency), the sample pointer rolls through the waveform at a rate of $Fs$ samples per second. To make the recording sound twice as fast, you can replay it at $2 \cdot Fs$. But, what happens when, to make the two waveforms catch up with each other, you have to run the "record" at $1.056 \cdot Fs$?

When the waveform is played back at normal speed, the sample sequence is $0, 1, 2, 3, 4, 5, 6, \ldots$. When the waveform is played back at twice its normal speed, the sample sequence is $0, 2, 4, 6, 8, 10, \ldots$.

But for non-integer factors, you would have to "land" in between two samples. For example $0, 1.056, 2.112, 3.168, \ldots$.

This is the reason for using interpolation. For example, to work out the sample value at time $1.056$ when you know the sample values at $1$ and $2$, you could use plain simple linear interpolation. This assumes, that when the Analog-to-Digital Converted (ADC) is not sampling the signal, this signal evolves smoothly (as a line) between its values. For high enough $Fs$, this assumption might be valid.

So, suppose that your signal is in $x[n]$ and you want to work out the signal value somehwere (let's call this $k$) in between $n .. n+1$. So, $n < k < n+1$. Then, $v = x[n] + (k-n) \cdot (x[n+1] - x[n])$.

This method is implemented in numpy 's interp() function or even MATLAB's interp1().

Because of this need of having to read in between samples, interpolation is used in such operations as reverberation, chorus, pitch shifting and others. This linear interpolation method will work in general for high enough $Fs$ but as the speed-up/slow-down factor becomes bigger and potentially, applied over a longer block of audio, the artifacts become more audible. This is what more complex methods of applying pitch shifting (such as the phase vocoder), are trying to address.

Hope this helps.


Feature Extraction of FFT for One Class SVM I'm looking for a good way of extracting features from the frequency domain of vibration data for a one-class support vector machine.

The image below shows an example from the dataset found in this link.
The dominant peak found at around 960 Hz corresponds to the "main frequency" of the system. E.g rotation*number of gears.

Im keeping this question fairly open ended, because i know there are several methods for doing this and I'm mainly looking for your suggestions and experience.

Would it be possible to analyze this so I get some specific information regarding these peaks caused by the faulty gear tooth?

Assuming you have epochs/segments of data. For this kind of signals, it is a safe approach to extract features using wavelet representations.

Using FFT might work as well, but I dont know how problematic would be the stationarity assumption in this kind of applications. Besides, FFT estimates for this kind of signals is sometimes very noisy. If you insist on this kind of representation, I would suggest estimating the power spectral density using Welch or multitaper methods. These usually deliver a nicer spectral estimation.

In any case, using any spectral features as direct input to your SVM might not be the best idea. Irrespectively of what features you extract (FFT, Wavelets, etc..) I would include a data dimensionality reduction step (PCA would be the first thing I'd try)


YIN pitch estimation algoritm-simplified explanation can someone please explain the steps of YIN pitch detection algorithm in a simple way "especially the last 3 steps".

here is the research paper of YIN algorithm
http://audition.ens.fr/adc/pdf/2002_JASA_YIN.pdf
Thanks in advance

Edit1(as i understand):
I understand the autocorrelation to find out if the signal periodic or not "step1", then we will take the difference between the original signal with zero lag"delay" and subtract it with a delayed one "step2",now we have the negative peaks but to make the algorithm works and pick up the first delayed negative peak with the lowest delay we have to remove the first negative peak with 0 delay"step3", the fourth step is to set a threshold to get the most negative peak(i am not sure if i am right)"step4",step 5 & step 6 i dont understand them,can you please tell me if i am explaining right or wrong

There is an older question about YIN where I am expressing lack of impressing.  The only thing that is novel in YIN is the Cumulative mean normalized difference
function (Step 3):

what is plotted above is the Average Squared Difference Function (ASDF) which is similar to the older AMDF.  the purpose of this Cumulative mean normalized difference function is to normalize the result (for the purpose of a threshold), and to do something to avoid picking the null at a lag of 0 samples from being picked as the lag for the best match.  they could have just said they won't start searching for the lag with the lowest null (or the null they want for the purpose of measuring the period) after the first lag that has a sign change for the slope.  that would be the first peak, and then they can start searching for the lowest null without worrying about picking the null at zero.

i'm not impressed with YIN.  i think i have a better idea for avoiding picking the wrong peak.

here's another one.

think of the Autocorrelation function as an upside-down ASDF.  so then you're looking for peaks instead of valleys.

UPDATE:

as per request, Step 5, parabolic interpolation of peaks (or valleys) is discussed in this answer and the Julius Smith reference therein.  it is for getting sub-sample precision from the ACF evaluated at discrete lags.

step 6 is about avoiding picking the wrong peak and avoiding indecisively switching back-and-forth between candidate peaks.  this is usually what we call "the Octave Problem".  i discuss it in this answer.  the YIN paper one identifies the problem and does not satisfactorally deal with it.


How are pulse trains generated in a radio transmitter? When simulating a PSK signal in MATLAB or Python, I often create an impulse train where each impulse is modulated with a symbol. This modulated impulse train then passes through a pulse shape filter to yield the PSK signal.

However I am wondering if this is indeed how PSK signals are generated in real radio electronics hardware. Is there an electronics module which generates modulated impulses which passes through an analog filter? If so, how narrow are such pulses typically?

I know the block diagram representation of PSK implementation in hardware often shows dataflow split into two paths, one for In-phase and one for Quadrature-Phase. Each path passes through a Root Raised Cosine filter before combining back again. However it is not clear to me how this RRC module is implemented in hardware.

Arbitrary Waveform Generators are a special type of test equipment that can generate, as the name implies, arbitrary waveforms as output signals that can be fed into transmitters or other circuits. They vary significantly in sophistication.


Convert frequency response to impulse response I am running an acoustic simulation.  I put in a frequency and get out a complex number; the frequency response for that frequency.

I wish to generate a 1024-sample impulse response (supposing my output sample-rate is 44.1kHz)

How do I go about it?

I think I have the first step: take the Nyquist frequency 22.05kHz, and cut it into 512 equal slices.  So each slice will have a bandwidth of 22050/512 Hz

Now I feed these frequencies into the simulator and get back 512 complex numbers; the first will be 0+0i corresponding to 0Hz

Now I'm following an instruction that says: 'consider 0Hz now to be the symmetry point, and add 511 negative frequencies.

So I have
    {
        Re(511) + i.Im(511),
        ...,
        Re(1) + i.Im(1),
        0 + 0i,
        Re(1) + i.Im(1),
        ...,
        Re(511) + i.Im(511),
    }

1023 frequencies.

Now I'm confused. I don't know what I'm doing any more.

Apparently I need to IFFT this, and that is going to give me an 1023 sample impulse response?

Could somebody help me get my feet back on the ground?

What would the MatLab code for this look like?

(Also, is there any way to do it from the Bash shell using Python or something similar?)

Here is a MATLAB example. Inefficient code but hopefully easy to read

%% Example: first order bandpass from 200 to 2000 Hz sampled at 44.1 kHz
n = 1024; % FFT length
fs = 44100; % sample rate
f0 = [200 2000]; % bandpass range
[b,a] = butter(1,2*f0/fs); % these are the difference equation coefficients
b0 = b(1); b1 = b(2); b2 = b(3);
a0 = a(1); a1 = a(2); a2 = a(3);  % for easier readability

% Step 1: build normalized frequency vector
om = 2*pi*(0:n/2)'/n;

% step 2: calculate the frequency response at every frequency
z = exp(j*om); % this is the "z" of the z-transform
% calculate the Z transform ,
H = (b0+ b1*z.^-1 + b2*z.^-2)./(a0+ a1*z.^-1 + a2*z.^-2);

% Step 3: make it conjugate symmetric
H = [H; conj(H(end-1:-1:2))];

% Step 4: inverse FFT
h = real(ifft(H));

% Step 5: verify against direct evaluation of the time domain difference equation
d0 = zeros(n,1); d0(1) = 1; % unit impulse
href = filter(b,a,d0);
err = 10*log10( sum((h-href).^2)./sum(href.^2));
fprintf('Error = %6.2f dB\n',err);


Circular Convolution and FFT of power 2 I am currently trying to understand the Chirp Z Transform, and compute an FFT for a vector that has size which is not necessarily a power of 2. The key step in the Chirp Z transform to recognize that the boils down to computing a circular convolution of two vectors of length $n$.
My question is: How does one compute the circular convolution of two vectors say ($x[1],x[2],...,x[n]$) and ($y[1],y[2],...,y[n]$) using DFTs only sizes that are a power of 2? I know that by the convolution theorem, I can find the circular convolutions of two vectors applying the DFT (and it's inverse) a few times, but I only have power of 2 sized DFTs at my disposal to compute this for arbitrary $n$. I considered padding the vectors to increase their size to the nearest power of 2, but it seems as if one must carefully construct such a padding to get the right result (since we need to compute modulo $n$ and not $2^k$).
For more context: I was going through Dilip Sarwate's answer on this topic at this link: https://math.stackexchange.com/questions/104148/chirp-transform-and-convolution[on math stackexchange]1 and was unable to figure out how one can compute a circular convolution efficiently.

Circular convolution is just linear convolution aliased by DFT length $n$. The length of linear convolution of $a$ and $b$ will be $2n-1$. So take $FFTs$ of $a$ and $b$ , padding each of them to length nearest power of 2 more than or equal to $2n-1$. Multiply the corresponding $FFTs$ point by point to get a power of 2 length sequence and take $IFFT$ of it. This sequence is actually the linear convolution of $a$ and $b$ since we had done enough padding before taking their individual $FFT$. Let this sequence be named $c$. Now, alias in time domain by shifting copies of $c$ by $n$ and adding them on top of each other.
$$
d[m] = \sum_{k=-\infty}^{k=+\infty} c[m - nk]
$$
The final output you want is $d[m]\,\text{for}\, 0\le m \le n-1$


Remove multiplied white noise from periodic signal I need to remove noise that is multiplied to a periodic signal using a maximum of 7 periods having no information about the noise.
I have tried to use auto-correlation:

 abs(ifft(fft(signal.*noise).*conj(fft(signal.*noise)))),

But I'm losing too much information from the original signal.
Any suggestion to recover the original signal keeping maximum points index information?

If it helps, the data set I'm using is here (Matlab mat file)

I suppose you can use synchronous averaging in this case due to the signal being periodic and the nature of the noise (as long as it is sufficiently uncorrelated). See MATLAB example:

ts = 1e-4; % sampling time
tend = 10; % end time
t = 0:ts:tend-ts; % time vector

fx = 10; % periodic signal frequency
x = 5*sin(2*pi*fx*t); % periodic signal

switch 2
    case 1
        n = 1*rand(1,length(x)); % noise in (0,1)
    case 2
        R = 50; % upsampling factor
        nn = 1*rand(1,length(x)/R+1);
        tn = 0:ts*R:tend;
        nnn = timeseries(nn,tn);
        nnnn = resample(nnn, t, 'zoh');
        n = nnnn.Data(:).'; % noise that looks like OP's plot
end

%y = x + n; % additive
y = x.*n; % multiplicative

nsp = 1/(fx*ts); % no. samples per period
NP = fx*tend; % no. periods

avg = zeros(1,nsp); % initialise
for k = 1:NP % loop over number of periods available

    Is = (k-1)*nsp+1; % index to start of period
    Ie = k*nsp; % index to end of period

    avg = avg + 1/NP*y(Is:Ie); % add period to average
end
avg = 2.*avg; % scale dependent on noise and how it enters the signal

plot(t(1:nsp),y(1:nsp),t(1:nsp),avg,t(1:nsp),x(1:nsp)) % plot signals
legend('signal with noise', 'average', 'signal')


Bluestein's algorithm to evaluate the DFT from $f_o$ to $f_o + k\Delta_F$ Briefly, the convolution between $x(nT) e^{-j2\pi f_o nT} e^{-j \pi \Delta_F Tn^2}$ and $c(nT) = e^{j \pi \Delta_F T n^2}$ multiplied $e^{j \pi \Delta_F T k^2}$ allows me to find the DFT $X(f_k = f_o + k\Delta_F)$.
Let $N_x$ be the number of samples of $x(nT)$, I found that $c(nT)$ must be at least of length $2N_x$ but I don't know why. Can someone tell me how do I choose the number of samples of $c(nT)$?

The convolution in Bluestein's algorithm is a linear convolution, not circular. Since the approach is to use standard FFT's for the convolution operation, the underlying waveforms being convolved must be zero padded out to $M \ge 2N-1$ to compute the linear convolution accurately.


OFDM modulator in LTE My understanding of what an OFDM modulator does, is that it takes a parallel sequence of QAM symbols of size M, pads it with zeroes corresponding to the zero subcarriers, then passes them through an IFFT block of same size M. IFFT is used because it can digitally create subcarriers.

My questions

1- Is the result output of the IFFT operation, an OFDM symbol of length M too?

2- Is EACH OFDM symbol going to be weighted sum of the ALL the QAM symbol?

3- How do things change if we have two antennas instead of one, do we have an OFDM modulator to each antenna? For example when using the Alamouti scheme? I have no clear vision how things change in this case, i.e what is the input at the IFFT now?

I am looking forward for answers thank you.

Definitely yes. But after IFFT operation data block is CP extended before transmission.
Each sample of OFDM symbol from the output of IFFT is weighted sum of all QAM/PSK modulated subcarriers from the input of IFFT.
In the MIMO case all spatial-time algorithms are applied in the frequency domain. So in transmitter you should do encoding with QAM subcarriers before IFFT and do the distinct IFFT for each spatial channel created. While in receiver you should do decoding after FFT applied separately to each spatial channel again. By doing this you exploit the fact OFDM is multiplexing of the large number (actually M) of narrow band channels and by appropriate processing (actually FFT/IFFT pair and CP insertion technique) we can treat each subcarrier as narrow band. It is very useful in the case of MIMO because MIMO decoding algorithms are much too simpler when dealing with narrow band (so flat fading channel is the case) and for each subcarrier we can define simple MIMO matrix of weights $[W_{ij}]$. That's why OFDM + MIMO is so popular combination today.


How to do a Fourier transform on input from a microphone and obtain the average amplitude of a certain frequency range I'm building a sound visualizer with LED strips controlled by python, and I'm trying to learn how to use pyAudio to perform rapid Fourier transforms on live mic input and obtain the average amplitude of a certain range of frequencies.

How do I open a sound stream and consistently update the current amplitude (say, a variable called "amplitude") of a given frequency range with the live FFT data? I have set up all the hardware and all the basic pyAudio functions work so far.

If you know the frequencies you're interested in, you might be better off just calculating the power / energy at those frequencies rather than using a full FFT.

The best algorithm to do that is probably the Goertzel algorithm.  There is a python implementation of that here, and a related SP.SE question here.


Transport Block calculation in LTE Downlink I would like to understand the relation between the TB size calculation in the TS36.213 document and in the simulator (https://www.nt.tuwien.ac.at/research/mobile-communications/vienna-lte-a-simulators/). I’m going to make an example so maybe you can help me to see this relation:

Method 1 (from the TS36.213)

Calculation Procedure for downlink(PDSCH) is as follows :

i) refer to TS36.213 Table 7.1.7.1-1

ii) get I_TBS for using MCS value (Let’s assume MCS is 1. in this
case, I_TBS is 1 )

iii) refer to TS36.213 Table7.1.7.2.1

iv) go to column header indicating the number of RB (Let’s assume
that RB is 50)

v) go to row header ‘1’ which is I_TBS

vi) we would get 1800 (if the number of RB is 50 and I_TBS is 9)

vii) (This is Transport Block Size per 1 ms for one Antenna)

And this is method 2 (from the Vienna LTE System Level simulator):

i) This is the formula:

TB_size_bits = max(8*round(1/8*(the_RB_grid.sym_per_RB_nosync .* num_assigned_RB .* modulation_order .* coding_rate * 2))-24,0);
ii) if the_RB_grid.sym_per_RB_nosyn= 80

iii) num_assigned_RB=50

iv) modulation_order = 2

v) coding_rate=0.0762

vi) That gives TB_size_bits = 1195 bits

Do you know how I can go from one method to another or if there is any relation?

Thank you. Best regards,

Natalia.

I found this ppt: http://www.ece.drexel.edu/walsh/Gwanmo-Nov11-2.pdf (slides 9-12) where I saw a bit of light regarding on how the method 1 calculates TB size:

For data, there are:

10 OFDM symbols per subframe x 12 subcarriers = 120 (we are counting the 6 RS symbols and excluding the PDCCH symbols)

Therefore,

120 symbols x 50 RBs x 2 bits/symbol (QPSK) = 12000 bits

For QPSK we know that the efficiency is 0.15 (approx.), so:

12000 x 0.15 = 1800 bits, as in table Table7.1.7.2.1 3GPP TS 36.312 for 50 RBs and I_TBS=1 (sorry, in my question I made a mistake and wrote I_TBS 9 the second time: can’t edit it in the forum)

However, in the simulator, using method 2, we are considering:

Number of non-sync symbols per RB= 12×7 – 4 (RS symbols)=80 OFDM symbols per slot

80 non-sync symbols x 50 RBs x 2 bits/symbol x 0.0762 x 2 slots =1219 bits

1219 bits -24 bits for the CRC = 1195 bits :)

Please, let me know if you don’t agree with these methods.

Thanks,

BTS.


What are some current algorithm(s) for finding if two images are copies of each other? Could anyone enlighten me on what are the current best ways to detect if two images are the same (i.e. a copy or edited copy of each other)? What feature detection method do they use (SIFT, DWT, pixel intensity, etc.)?

I've read a lot of algorithms that focus on detecting similar objects and so on, but I am interested in reading on some for detecting whether two images are copies of each other. I've done some googling and most algorithms/papers I've found are pre-2014 or are about detecting similar images, not exact copies.

This is not a simple problem. You need to define these 'simple edits' and check for all of those.

For example if you want to find out brightness changes, you cannot effectively use template matching, you would have to go with feature matching (contour, keypoint etc.).

For a general solution, you need to extract a 'signature' from an image which for instance can be done with a combination of histogram and feature points (edge, corner, blob, Shi-Tomasi, LoG etc.)


Inverse Hilbert Transform The Hilbert Transform of a 1D/real-valued vector signal returns the analytic signal, x, from a real data sequence, xr. The analytic signal x = xr + jxi has a real part, xr, which is the original data, and an imaginary part, xi, which contains the Hilbert transform.
hilbert uses a four-step algorithm:

1. Calculate the FFT of the input sequence, storing the result in a vector x.

2. Create a vector h whose elements h(i) have the values:

     1 for i = 1, (n/2)+1

     2 for i = 2, 3, ... , (n/2)

     0 for i = (n/2)+2, ... , n

3. Calculate the element-wise product of x and h.

4. Calculate the inverse FFT of the sequence obtained in step 3 and returns the first n elements of the result.

This algorithm was first introduced in [8]. A python implementation of this can be seen bellow:
from scipy import linalg, fft as sp_fft
import numpy as np

def hilbert(x, N=None, axis=-1):

        x = np.asarray(x)
        if np.iscomplexobj(x):
            raise ValueError("x must be real.")
        if N is None:
            N = x.shape[axis]
        if N <= 0:
            raise ValueError("N must be positive.")
        print(x.shape,N,axis)

        Xf = sp_fft.fft(x, N, axis=axis)
        print(Xf.shape)
        #plt.plot(Xf)
        #plt.show()
        h = np.zeros(N)
        #plt.plot(h)
        #plt.show()
        if N % 2 == 0:
            h[0] = h[N // 2] = 1
            h[1:N // 2] = 2
        else:
            h[0] = 1
            h[1:(N + 1) // 2] = 2
        print(h)

        if x.ndim > 1:
            ind = [np.newaxis] * x.ndim
            ind[axis] = slice(None)
            h = h[tuple(ind)]
        x = sp_fft.ifft(Xf * h, axis=axis)
        return x

This code was taken form Scipy implementation Scipy.signal.hilbert. I am looking to invert/reverse this process (inverse_hilbert), the best description i have found to do this is from Mathworks Inverse Hilbert Transform

However the real and complex arrays currently have me at a loss not sure how this feeds into this equation, or if this is the correct equation as wikipedia has a different equation to my understanding.

If we create a random complex signal and compute the hilbert Transform on it we get the following 2 arrays one real and one imagery. I am looking to reverse this transform any help would be appreciated.
import matplotlib.pyplot as plt

x = np.arange(0, 30, 0.1);
y = np.sin(0.05*x)+np.sin(6*x)+np.cos(3*x)
plt.plot(x,y)
plt.title('Complex Waveform')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.show()

x_a = hilbert(y)

plt.plot(x,x_a.real, label='Hilbert Real', alpha=0.5, lw=2)
plt.plot(x,x_a.imag, label='Hilbert Imag', alpha=0.5, lw=2)
plt.grid()

Any comments here to help my under standing would be appreciated, really looking for a step by step way to reverse this transformation.
[8]: Marple, S. L. “Computing the Discrete-Time Analytic Signal via FFT.” IEEE® Transactions on Signal Processing. Vol. 47, 1999, pp. 2600–2603.

I think the confusion comes from the fact that the command hilbert in Scipy (and also in Matlab/Octave) does not just compute the Hilbert transform, but its output is the analytic signal. So if $x(t)$ is the (real-valued) input to such a function, its (complex-valued) output is
$$y(t)=x(t)+j\mathcal{H}\{x(t)\}\tag{1}$$
Clearly, if you want to obtain $x(t)$ from $y(t)$, you just need to take its real part.
Concerning the inverse Hilbert transform, as you've already found out, the following holds:
$$\mathcal{H}^{-1}\{x(t)\}=-\mathcal{H}\{x(t)\}\tag{2}$$
The relation $(2)$ holds no matter if $x(t)$ is real-valued or complex-valued. However, if you use the function hilbert on a complex-valued signal, you throw away information that cannot be retrieved, because you throw away the negative frequencies, which are not redundant for complex-valued signals (unlike for real-valued signals). So the function hilbert  cannot be inverted if the signal is complex-valued.


Historical notes of filter designs? Anyone know what references might have information about the history of any particular filter design?

Since seeing just the equation or implementation hinders the process by which the filter (such as the biquad filter) was actually discovered.

Along the line as Fat32's suggestion to consult IEEE papers, you might try to obtain a copy of "Digital Signal Processing" (IEEE Press Selected Reprint Series), Edited by: Lawrence R. Rabiner and Charles M. Rader, published by: IEEE press (1972).

That book contains reprints of many early (starting around the middle 1960s), and influential, scholarly papers on digital filtering. And of course those early papers contain references to even earlier papers.


Effect of Hanning Window on phase of FFT I'm calculating the beat of song through the frequency domain. I have 60 seconds of audio sampled at 48000Hz, multiplied by a Hanning window in the time domain. This gives me a resolution in the frequency domain of 1BPM, and I'm interested in frequencies from ~10 BPM to ~300BPM, which is only the first 300 bins of the 48000*60 bin length FFT. How does the Hanning window affect the phase of the song's beat? Using no window (really a rectangular window) results in a different phase than the Hanning window does.

It't really more of a question of what the rectangular window is doing to your phase relationships, because it doesn't deemphasize the "edge effects" that occur at the beginning and end of your block of samples, which is the whole point of using a non-rectangular window.

Try this experiment: Do a number of FFTs with lengths ranging from 59*48000 samples to 60*48000 samples, using both window types. You'll find that the Hanning window gives much more consistent results than the rectangular window does. The two results will agree best when the length is close to an integer number of beats.


How to choose one among different types of LPF? My signal graph is as follows

I have < 10 reputation, thus cannot post the photo...

The graph are simply roughly the series of -5, 0, +5, 0, -5, 0, +5, 0.... but with sawtooth-like small peaks around the prominent peaks.

I am going to do the peak detection. But what I need to count is the "prominent peaks".

i.e. The "sawtooth" small peaks around the main peaks should not be counted.

I have known that the small peaks are due to high-freq noises.

My range of interest is 0~4Hz.

Thus, before counting the peak number, I have to low pass filter the signal first.

Being a beginner in signal processing, I have no idea which type of filter is the most proper. Chebyshev or Bessel or the classic Butterworth?

Thanks in advance!

There are 4 main types of LP filters, Butterworth, Chebyshev1, Chebyshev2, and Elliptic filter. Each differs from the other in the oscillations it has, and in which part of the spectre it is.

Butter has the worst amplitude cutoff, eg the widest transition region, but no osclillations;
Cheb1 has better transition range, but has oscillations in the pass region;
Cheb2 is the opposite of Cheb1, has oscillations in the reject region, but not in the pass region;
Elliptic has the narrowest transition, so therefore the best, but has oscillations in both the pass and no-pass zone.

Heres a picture to help you visualize:

http://i.stack.imgur.com/rJMQw.png

Now, seeing as you need a peek detector, you should use a flat filter for the pass zone, so Butter or Cheb2. I advise Cheb2, because you don't really care about the higher frequencies.

In matlab, the function is cheb2ap, and is used like this:

[z, p, k] = cheb2ap(N, La)
% where La is the allowed ripple in the stopband (usually around 0-1.5 dB).

% And now the filter transfer function is acquired with:
b = k*poly(z);
a = poly(p); % b is the polynomial for zeroes, a is polynomial for poles
[H w] = freqz(b,a, N)  % N is the number of points for calculation

Now you have frequency response of your filter, H, and vector for its frequency, w, and you can plot the filter with plot(w, abs(H)).

Hope this helps :)


A general conceptual understanding of moving average I have a conceptual question here. I understood how to calculate a Moving Average(MA).

For example :

x(n) = 0, 1, 2, 3, 4, 5

MA = 0, 1/2, 3/2, 5/2,7/2, 9/2, where for calculation of the MA we only take the past values without any update

Is there also a case , where we use the newly calculated values and use that to update the MA for the next time instant?

x(n) = 0, 1, 2, 3, 4, 5

MA (n=0) = 0

MA (n=1) = (0 +1) /2 = 1/2

MA(n = 2) = (2 + 1/2) /2 = 5/4

MA(n = 3) = (3 + 5/4) = 17/4 and so on ?

Does this find an application anywhere?

Please help me with this understanding.

That system you're referring to is a first-order recursive filter

$$y[n]=\alpha x[n]+(1-\alpha)y[n-1],\qquad 0<\alpha<1\tag{1}$$

where $n$ is the (integer) time index, $y[n]$ is the output, and $x[n]$ is the input. You chose $\alpha=\frac12$.

That filter is also called exponentially weighted moving average. It has theoretically infinite memory and weighs the past input values with exponentially decreasing weights.


Digital filters in LabVIEW I am very new in signal processing and using digital filters. I have to use a low-pass filter to analyze my data in LabVIEW and have a question about it. Any help and advice is appreciated.

I am trying to simplify my problem here:

Let’s say there is a digital sine wave (made by LabVIEW) with $V_{offset}=1 \ \mathrm{V}$, $V_{peak}=0.1 \ \mathrm{V}$, $f=10 \ \mathrm{kHz}$, $N=2000$ (number of samples), and sampling rate $f_s=200 \ \mathrm{kHz}$. Now, if I pass this signal through a low-pass filter with cutoff frequency $f_c=1 \ \mathrm{kHz}$, then the output should be a constant number equals the DC offset (here $1 \ \mathrm{V}$), is it true?

Another question is the concept of “cutoff freq” and “sampling freq” as the inputs of the filters in LabVIEW. Cutoff frequency as an input of a filter makes sense to me but what is that “sampling freq” ? Can anyone explain to me please? I am very confused. Is it the same rate at which the sine wave is created? For example in the attached code, what is the real cutoff frequency (with $f_l=200000$ and $f_l=1000$)?

I have attached the screenshots of the Front panel and Block diagram of my simple vi.

Let me answer your two questions in turn:

For your first question, generally, yes that is correct; if you filter a 10KHz sinewave that has a DC offset with a filter that has a cutoff frequency below the frequency of the sinewave, then the sinewave would be rejected. The amount of rejection specifically depends on the performance of the filter, but given you said you have a 1KHz cutoff frequency, the sinewave is significantly higher and therefore sufficiently rejected.

I see in your plot that the order of the filter is 5, which for a Butterworth filter as also shown would have a rejection of 20dB/decade *5 (where 5 is the order of your filter), or 100 dB per decade. Depending on other factors such as your digital dynamic range, this suggests that you would be able to filter your 10KHz sine wave up to 100 dB (10KHz is a decade above the cutoff frequency).  The DC signal, which is below the cutoff frequency would pass through to the output, unless something in your system blocked DC or introduced other DC -offsets (which is possible).  Also the filter itself can have gain or loss, so the actual DC output level if it did pass through can be modified by this gain or loss accordingly.

For your second question, sampling frequency is the sampling rate for the signals passing through this digital filter implementation.  From the figure, you are using a sampling rate of 200KHz, and yes this would be the sampling rate of the sinewave that is created.

I hope this helped to clear up some of your questions.


Can I determine a system's $z$-domain transfer function from its pole-zero plot? Is it possible to generate the z-domain transfer function from a pole-zero plot diagram?

Without knowing specifics (ignoring a proportionality constant we'll call K as inidicated in the more compete answer above).  Note that K does not affect the dynamic behavior of the system, so it may not be of particular interest.

it goes like this...

H(z) is a rational function of the form
$$H(z) = \frac{(z-Z1)(z-Z2)...(z-Zm)}{(z-P1)(z-P2)...(z-Pn)}$$

where Z1,Z2...Zm are your zero values and P1,P2...Pn are your pole values.

I assume you know the exact values of you're poles and zeros. If you don't, you'll have to determine approximations from the PZP.


Hardware advice for machine vision project I am working in facial feature detection in video images. I was hoping to get some general advice about camera/lens technologies to determine if I can upgrade my camera/lens.

Basically, the problem is this:

I want to obtain high resolution images of facial features.

I currently use a Thorlabs DC1545M CMOS camera and a 16mm focal length lens. Generally, when I extract the facial features, the image regions have acceptable resolution if the face is about 40cm or less from the camera. If face moves further away from the camera, the resolution of the facial features starts to get too low for my purposes. Additionally, the camera/lens has a short depth of field so the face will often go a little ‘out of focus’ as the face moves.

For this region I was hoping to improve my camera and lens setup. However, I am not strong on camera hardware so I am not entirely sure if what I am thinking is available.

Basically, my ideal camera would do the following:

It would have an optical zoom.
The zoom would be programmable.
It would be able to have a long depth of field so the image stays in focus. Basically the face will be closer than about 100cm to the camera so I would like face to stay in focus.
I would like the camera to be sensitive to infra-red light (Near IR).

One other thing. Right now, I could use off-the-shelf camera/lens. However, I would also be interested in board level camera/lens components.

I was hoping someone could give me ideas about what I should be looking for.

Are there problems with my specification?
Do you know of any cameras/lens with the specification?
Do you have any general advice about what I should be thinking about before looking to buy a camera for this kind of research???

The things you can control with the camera selection are:

Resolution
Sensor size
Sensitivity
Frame rate

The things you can control with the lens selection are:

Focal length
Minimum depth-of-field (DoF)
Field of view (FoV)
Minimum aperture
Zoom functionality

Each of these interact with each other so it is often difficult to find the right trade off.

The sensor size and resolution are usually related. For example, it is hard to find sensors with resolution above 1024 * 1200 in a 1/3" sensor.

The resolution and thus sensor size will dictate what lenses can be used with the camera, then combined with the FoV of the lens will dictate how many pixels the face covers in the image.

First decide how far off-centre the face can move. That will guide you minimum FoV. The work out how far away the face can be, and how many pixels wide it must cover. Using the average width of a face, the FoV angle, and the depth you can work out the minimum resolution you need.

Using the minimum resolution and DoF, choose a camera and lens combination. Pay close attention to the maximum sensor size the lens can accomodate, and note that the FoV changes with the sensor size.

The DoF is a function of the lens focal length and the aperture. To increase the DoF you have to reduce the aperture. E.g. f/1.8 is wide open, f/22 is almost closed. A smaller aperture will let less light in, and thus the camera exposure and gain will need to be increased. Choosing a camera with high sensitivity will allow you to get better images with less gain increase while keeping the exposure short enough to not be blurry.

Motor driven zoom lenses are available.

Near-IR is possible with many camera manufacturers, as they actually have to put in a filter to block IR, so making a near-IR one is easy - just take the filter out.

But probably my best piece of advice is to contact the camera manufacturer and ask them. They have people that will take your specifications and recommend a camera / lens pair for you. Personally, I have used Point Grey cameras.


Separation of signals without loss of information about the signals themselves There is a signal that contains:

Rectangular pulses on frequency of 1 kHz
Rectangular pulses on frequency of 1.4 kHz
Other noises.

Required:

Separate 1kHz pulses from everything else and measure their amplitude and duration
Separate 1.4kHz pulses from everything else and measure their amplitude and duration

Was tried:

Use a positive IIR filter of the 6th order (at a frequency of 1kHz and at a frequency of 1.4kHz). After the filter, information about the initial amplitude and duration was completely lost.
Use comb filter. The output signal is smeared due to the mismatch of the input signal period and the filter period.
Searching and asking a question in my native language (Russian) did not help. I did not find any examples. The only possible solution proposed was the use of correlation analysis. But I could not work out, how to use it for solving my problem, again there was not a single explanation or example.

P.S. I’ve already asked this question, but maybe I wrote something that is not clear (I am not a native speaker of English), please ask, I will try to explain more clearly. I have one last hope on the international site that someone has faced this problem and will be able to share a solution or explain how to solve it. Thank you very much in advance!

P.S.S

  why do you want to filter? your signals look clean

In that question, the picture was just for example. In this question, the picture looks like my real signal.

  Yeah, classical XY problem: You're asking a strange question (you want to filter, but you want to have none of the effects of filtering?) without telling us why you want to do it, for what purpose, in what greater context.

In this thread, I tried to describe in more detail what I want.

  If I understand correctly and you want a system that picks some pulses and rejects others based on their fundamental period: a filter is not going to do it, at least not by itself.

I need to separate the 1 kHz pulses from everything else and measure their amplitude and duration, and then separate the 1.4 kHz pulses from everything else and measure their amplitude and duration.

  Does the amplitude stay constant for the duration of the pulse groups? Do you have a subset of possible pulse groups: number of pulses at the 1.4 KHz rate, number of pulses at the 1KHz rate? Please provide all info in terms of what these pulses are doing and why you are doing this may also help answer those questions you may not be sure how to ask. The sentence "at certain points in time it can change by 200 us..." was particularly confusing- what is "It", the duration? What does it mean for the duration to change by 200 us, increase to 300 us?

No, the amplitude can vary. I tried to describe everything in more detail in this topic.

  Do the pulses add (i.e., if a 1.4kHz pulse and a 1kHz pulse land on top of each other is the result the sum of the pulse heights)? Do you have to catch each and every pulse, or can you miss some?

Yes, when they fall on each other it will be the sum of the heights. The signal above shows this. In the ideal case, I would like to catch every impulse, but in fact you can also skip.

This is probably going to take multiple steps.  As a start, here is what I'd do:

You want a filter that retains pulse information, but filters noise.  The most extreme you can do with this is a rectangular filter that's $100\mu\mathrm{s}$ long, or possibly a raised-cosine filter that's about the same length.
Assuming that the pulse timing is very consistent, phase lock to each pulse train.  This means you have one 1kHz PLL, and one 1.4kHz PLL.  When you get ready for this step you may want to ask a separate question, but if your pulses always start with the same timing (i.e., if a "long" pulse just ends $100\mu\mathrm{s}$ later than "normal") then you can lock onto the leading edge of each pulse.
Any time the PLL phases indicate a pulse will overlap, ignore both -- basically, output "null" on both channels.  You probably need to allow room for those $200\mu\mathrm{s}$ pulses you speak of.
Detect long pulses by looking for unexpected energy $100\mu\mathrm{s}$ later than the expected center of the filter output.
Detect pulse height by looking at the height of the filtered pulse where the PLL tells you.

The very beginning part of development for this would be to record your data into a file so you can play with it, run it through your best guess at a matched filter, and look at it.


Obtain transfer function from signal flow graph I have a pretty complex signal flow graph with all the blocks and the branches completely describing an analog filter, in the Laplace complex variable s. The signal $X(s)$ is the input and $Y(s)$ is the output.

Is there any (free) tool which can analyze such a graph and provide the overall transfer function $H(s) = Y(s)/X(s)$?

Googling I found this Matlab tool, but the tool it is not free; anyway is there anything similar, in Matlab or in other environments?

I am afraid that for complex systems finding such free standalone tool would be hard. Obviously since you mentioned MATLAB, you should use Simulink then or this tool:

Fully Automated Flowgraph Analysis Tool. Manual to be found here.

Alternatively you might look into this Mason.m tool. Or this masonry tool.

I know that these guys also did similar work once.


Problem understanding the Expectation Operator I know that the Expectation Operator $E\{x\}$ four discrete values is
$$ \sum_k \alpha Pr(x = \alpha_k)$$

and its very intuitive when speaking out a formula which contains the Expectation Operator. But I often have some troubles when trying to apply it on real examples with numbers.

E.g. the autocorrelation is completely obvious when given as sum like:
$$ r_x(\tau)=\sum_n x_n x_{n-\tau}.$$ However, often principles like autocorrelation, correlation, cross-correlation etc. are formulated with the $E$ Operator. For the autocorrelation:
$$ r_x(k,l) = E\{x(k)x^*(l)\}.$$
I assume this can be written as:
$$ r_x(\tau) = E\{x(k)x^*(k-\tau)\}.$$

The product is the same as in the sum above. However, the $E$ Operator implies that I need to know some probabilities. How can one quickly see that this is the same as above with $\alpha$ is somehow $x_n x_{n-\tau}$ and $Pr=1$? This can't always be the case. Otherwise one could always write the sum instead of the E Operator, couldn't he?

I think about the E Operator as if it is needed for different random processes. This would imply that I can't formulate e.g. a cross-correlation as a sum over all products between two signals, as I used to do. A formula like $r_{xy}(k,l)=E\{x(k)y^*(l)\}$ would then be a more general case and could not be applied on a simple numerical example, without making assumptions about the processes statistical properties.

Consider a discrete-time, discrete valued random process $X[n,s)$ whose samples $x[n]$ at integer index $n$ belong to range space $R_{X_n} = \{x_1,x_2,...,x_n\}$ of discrete valued random variables $X_n$ indexed by $n$, such that at each $n$, the value $x[n]$ is the result of the mapping by $X_n$ into $R_{X_n}$; i.e., $X[n,s) = X_n(s) = x[n]$, where $s$ is the sample output of the experiment performed at time $n$.

The expectetation of a particular random variable $X_n$ of the random process is
$$\mathcal{E} \{ X_n \} = \sum_{x_k \in R_{X_n}} x_k P_{X_n}(x_k) = \int_{-\infty}^{\infty} x f_{X_n}(x) dx $$
where those $x_k$ belong to the range space of the random variable $X_n$ and $P_{X_n}(x_k)$ are the probabilities of those values. Note that the integral definition is also available for those who would still like to use a continuous random variable approach with impulsive PDFs; $$f_{X_n}(x) = \sum_k P_{X_n}(x_k) \delta(x-x_k) $$

The correlation between two (real) random variables $X_n$ and $X_m$ is defined as
$$\mathcal{E} \{ X_n X_m \} = \sum_{k} \sum_{l} x_k x_l P_{X_n,X_m}(x_k,x_l) $$
where $P_{X_n,X_m}(x_k,x_l)$ is the joint probability mass function of the two random variables, and signifies the probability of the pair $(x_k,x_l)$ in the range space of the joint mapping.

When those two random variables belong to the same random process, then the correlation is defined as auto-correlation of the random process.

Note that the second sum in your post does not represent an expectation but rather a possibly ergodic estimate of it for a WSS random process.


Averaging Multiple Magnitude Spectra I'm currently trying to implement some of the methods found in this paper on intelligent equalization - http://www.aes.org/e-lib/browse.cfm?elib=16792 -

The first part of the process is to build a "target equalization curve" via the means of "averaging the spectra of all songs in the dataset".

In MATLAB, if I do the following.

mag_1 = abs(fft(file));
mag_2 = abs(fft(file_2));

The two magnitude spectra have different resolution and number of frequency bins right?  Does it make sense to just mean(mag_1, mag_2)?  I can't seem to find anything that goes over averaging multiple magnitude spectra and it isn't defined in the paper how they achieved the averaged magnitude spectra for the ideal target curve.

Any help would be appreciated.

If you're averaging multiple (two, in your case) spectral magnitude sequences your result will be valid, correct, as long as the sequences are equal in length.  Now whether the "averaged spectral magnitude" result gives you useful information, that you will have to decide on your own.


Opposite of wavelet transform? Wavelet transform gives good time resolution for high-frequency events and good frequency resolution for low-frequency events.

=> I want to have complete opposite of wavelet transform, where I get good time resolution for low frequency events? Is there any known transform or specific window to get it?

wavelet transform also serves to you what you want to do, all you have to do is to first apply a low-pass your signal so as to keep the frequency range you want to scan only, another approach is to transform Gabor where you have to define the size and shape of the analysis window, but I would recommend better  using wavelet transform applying a lowpass filter your signal first.


FIR vs IIR complexity IIR filters have feedback whereas FIR filters don't use feedback.
Then why is it said that FIR filters sometimes have the disadvantage that they require more memory and/or calculation to achieve a given filter response characteristic?

Can you please explain explain what is it that makes FIR filter computationally more complex compared to IIR although it doesnt use feedback?

while it is true that FIR filters don't necessarily use feedback, there are FIR filters that do use feedback.  the name for these kind of FIRs coined by Julius Smith and Avery Wang is "Truncated IIR" or "TIIR".  a moving average filter, where the moving sum is computed with an accumulator where the new sample is added and some old sample is subtracted, that is an example of a TIIR.

another method of doing an FIR is the use of the FFT and "fast convolution".  but there is no feedback with this method.

the measure of "complexity" of filters is usually how many multiply-accumulate (or "MAC") instructions, per sample, are needed to do the job.  if sharp filtering is needed, impulse responses must be long and, for a simple FIR, that means many MAC instructions.  but a low-order IIR can still have a long impulse response, even if there are few MAC instructions because of the feedback.

so if you can do the same basic job (like low-pass filtering) with a low-order IIR, in comparison to doing the essentially the same job with an FIR that would have scores or hundreds of MAC instructions, it should be clear what advantages regarding complexity the IIR method will have.


Entropy/frequency of a natural image I have a C# desktop application.

I am using a wrapper to OpenCV which is EMGU.

I am just performing a basic motion detection operation.

I am using component labelling and 2 frames differencing to determine motion changes.  I am also using an averaging background approach to update the background which i use to 'subtract' from the current and previous frames.

All this works well.

Now consider this image:

Because this bush is close to the camera, movement will be detected by my current methods when the wind blows.

Also, consider this image:

Here we see a cobweb to the right of the image. Again when the wind blows I get motion detection.

My possible solutions:

Look for shapes. In the case of elongated shapes like the cobweb where the height is greater than the width (to a certain scale) then ignore that movement.
But, I found I am missing motion when someone walks on the pavement at the back of the image from across the road.

Erosion.
Again I am missing genuine motion like people and bikes.

Compare Histograms
I have found it is not reliable.  A histogram will only measure the number of different colours (in its own bin).  It is possible for the colours to change due to the saturation and lighting or/and the same numbers of colours are found because another object has passed by.

So, I thought about using atrophy. I did post an earlier question to this here:

Entropy of an image

but as the person who so kindly explained entropy suggested that I should post as a specific/separate question here...

I want to explore whether the frequency or/and entropy would help me here?  Should I look at spatial or frequency distribution?  Would frequency be best to look into if I want to avoid colour changes?

Would either approach be faster than the other as I want it real-time as possible.

Any eduction for me here would be great.

Entropy would be a good metric in order to separate trees. However I would suggest you to calculate the entropy of gray values within a rectangular window. Consider this image:

After entropy filter, I could obtain this one:

Note that the trees (and also other informative areas) look white and it's possible to segment them out with a simple threshold. Try and see whether it also works for your case.

Due to the possible pre-computations, this could be optimized to run in real-time.


How can I connect multiple microphones to my computer (more than 4) and process the outputs independently? I am trying to create something like a heat map for sounds originating in a classroom to study the discussions inside the classroom. For this I am instaling directional microphones in various parts of the classroom and inputing the amplitude of sounds captured by them. I need to take the inputs from the microphones and process all of them separately. What hardware/software should I use for this?

PS: I have a mixer board but it'll give me a single channel input. I want separate inputs from all the microphones

Depending on what kind of application is that and what kind of microphones (dynamic, condenser, home-made electret ones) solutions can be different. Are you going to process everything in real-time or to make your analysis offline?

If your microphones doesn't require Phantom power and can be plugged into mic-in, then you can buy few of cheapest USB sound cards, plug them to your laptop, run your linux, create virtual device with JACK and you can record. But do not expect to run in real time - most likely there will occur lot's of problems with synchronisation (had them on my own).  So it is good for not a real-time application. Here is some article about that: Using Multiple Devices with JACK and topic from Audacity guys.
If you want to use microphones that require a phantom power, then you should buy/borrow some 4+ channel USB/FireWire audio interface. You plug it into your computer and record whatever you want. No problems with latency, everything is synchronised. What program can you use? That's up to you - on linux I suggest Ardour on Windows I am using REAPER a lot. But for scientific applications just use the ASIO drivers and record tracks.


What is the Uncertainty or error due to sampling frequency? If I want to detect the distance between two peaks in a signal sampled at 100 Hz. Then is it right to say that the first and second peaks occurred within 0.01 seconds of uncertainty? Therefore, the distance between peaks is 0.02 seconds uncertainty?
Then, based on that, how does the uncertainty propagates if I have a group of peaks and I want to find say the uncertainty in the mean and standard deviation of distance?
Does anyone have any idea what this uncertainty is called?

Uncertainty IS NOT determined by the sampling frequency (i.e., 100Hz) only.
Uncertainty is mostly determined by both the Bandwidth of the signal you sampled (as stated by Marcus Muller) and the sampling frequency.
In fact, the uncertainty between sampled peaks can be solved by properly designed interpolation (or DAC) as long as the sampling theorem is met. The main reason is depicted as follows.
Although Nyquist-Shannon sampling theorem state that no theoretical reconstruction loss as long as the sampling frequency is twice larger than the bandwidth. But, YOU NEED TO RECONSTRUCT THE SAMPLED SIGNAL before the temporal analysis. Most of the analyzers forget to reconstruct (interpolate) their sampled signal while analyzing.
This sample code sample a simulated 2 Hz analog signal,x_a, (with sampling rate 1500 Hz) to a discrete-time signal,x_d, with 5 Hz sampling rate (nearly Nyquist rate) and then reconstruct this back,x_r.
%% This code simulate the AD/DA processing discussed in Chapter 4.8.3 [1]
close all; clear all;

%% parameters.
% analog
analog_fps = 1500;
analog_window_time = 3; %sec
t = 0:1/analog_fps: analog_window_time-1/analog_fps;

% digital
digital_fps = 5;
n = downsample(t,analog_fps/digital_fps);

% ADC: Quantizer
X_m = 1; % Range
B = 10;% Bit number.

%% Signal generation
freq_hz = 1; % Hz.
x_a_1 = 0.5*cos(2*pi*freq_hz*t+0.1);

% add a small high-frequency component as asked.
signal_freq = 2; %Hz
x_a_2 = 0.5*cos(2*pi*signal_freq*t+pi/2);
x_a = x_a_1 + x_a_2;

%% ADC
% Sampling
x_s = downsample(x_a,analog_fps/digital_fps);

% Quantizing (abs of input value should not over 1)
% x_d = Quantizing(x_s,B,X_m);  % A For complete ADC, a quantizing should
% be added here.
x_d = x_s; % For basic case, we skip the quantizing here.

%% DAC
% up sample / DAC
x_up = upsample(x_d,analog_fps/digital_fps);

% LPF (Reconstruction Filters)
h = intfilt(analog_fps/digital_fps,4,0.9);
%% Important
% please not the parameter 0.9, ideally should be 1 for Nyquist rate.
% 0.9 here is the ratio of Nyquist.
% Given known limit band signal, a shorter ratio can enhance SNR by oversampling.
%  (i,e, here I filterout the freq larger than 2.5(Nyquist rate) * 0.9 = 2.25Hz)

x_r = filter(h,1,x_up);
x_r(1:floor(mean(grpdelay(h)))) = [];
x_r = [x_r zeros(1,floor(mean(grpdelay(h))))];

%% Display

figure;
plot(t,x_a);
hold on;
plot(n,x_d);
plot(t,x_r);
title('analog signal (1500Hz) v.s. digital signal (5Hz) v.s. Reconstructed signal (1500Hz)');
legend('x_a','x_d','x_r');

If you forget to interpolate, the uncertainty of x_d is very high!! However, as long as you properly interpolate the x_d, you will obtain the x_r where the uncertainty is much much smaller.
To sum up, as long as the band-limited signal is sampled higher than the Nyquist frequency, properly reconstruction (interpolation) can theoretically reduce the uncertainty to a very small amount (which depends on the interpolation filters you designed). The interpolation filter sample I provided above is FIR and zeros phase filter, which means finite impulse response (i.e., no error propogation and phase distortion).
DO NOT FORGET to RECONSTRUCT YOUR SIGNAL, or you will overestimate the uncertainty.


Beamforming with directional sensors Having photo-diodes with a very narrow field of view (+-10 degrees for half power) what should I consider when designing a beam-forming algorithm?

I know that I could use many photo-diodes mounted in an arc to cover 180 degrees for example but the uncertainty of 10 degrees is too big. As an alternative to classical beam-forming, what would be the methods to use to isolate a signal coming from one direction with a high precision (under 1 degree)?

My goal is to analyze the flickering of lamps, knowing their position from analyzing camera data, I want to isolate the signal coming from their direction.

Given the horizontal and vertical component of the distance vector from the source to the sensor array be X and Y, then the angle between the sensor array and the source is given by:
$$ \theta = \tan^{-1} (Y/X) $$

If the linear distance between each sensor is $d$, each sensor should experience a delay:
$$ \delta t = d \sin(\theta)/c $$

Now if you can measure the delays (high precision measurement will lead to higher accuracy), and have the knowledge of $d$, you can find the angle $\theta$.

Warning: Practical errors and considerations are not taken into account in this approach.


A question about Fourier transformation Hello, this is my first time actually asking in stackexchange. I am a computer engineering student and currently i am doing a linear system course (i don't really know how this is equivalent in engineering course labeling as my major is apparently in faculty of computer science but it's related a bit toward signal processing), i hope i wasn't in the wrong category

So I was to find the fourier transformation of the graph in the picture, and i need to plot the graph in matlab as well. I was allowed to look into the table but I couldn't find it hence I used brute force(integrals) to find the $X(\omega)$ of the signal, which i ended up finding but still there's $e^{-\omega.j}$ stuck which confuses me as i believe the plot should be 2 dimensional only (toward x and $\omega$, not toward j). can somebody help finding any faults/things i missed here?

Regarding your integral, it is not quite correct I think. You used partial integration to solve it, which works. Let's consider $\int t {\rm e}^{-\jmath \omega t}{\rm d}t$ and call $f(t) = t$ and $g'(t) = {\rm e}^{-\jmath \omega t}{\rm d}t$ so that $g(t) = -\frac{1}{\jmath \omega}{\rm e}^{-\jmath \omega t}$. Then using $$\int f(t) g'(t){\rm d}t = f(t)g(t)-\int f'(t) g(t) {\rm d}t,$$ we obtain $$\int_0^1 t {\rm e}^{-\jmath \omega t}{\rm d}t = \left[-\frac{t}{\jmath \omega}{\rm e}^{-\jmath \omega t}\right]_0^1 + \frac{1}{\jmath \omega}\int_0^1{\rm e}^{-\jmath \omega t} {\rm d}t.$$ It's close to what you did, but you used $g'(t)$ in the first term instead of $g(t)$.

  ended up finding but still there's e(−ω∗j) stuck which confuses me as i believe the plot should be 2 dimensional only (toward x and ω, not toward j).

Note that $j$ is not a variable! It's the imaginary unit $\sqrt{-1}$. Hence, your expression only depends on one variable, which is $\omega$. Matlab will recognize $j$ as the imaginary unit by default, so you should be all set for plotting it.


Uncorrelated but overlapping spectrum. This answer may be straight forward but I cannot figure out. One can understand the difference between cross-correlation and convolution from the link below:

What is the difference between convolution and cross-correlation?

The question is can there be two complex signal which are uncorrelated but their spectrum overlaps? In other words spectrum is not disjoint or product of power spectrum of each signal is not zero in each frequencies?

If you mean correlation at the origin (no signal shifts), then yes so long as all the overlapping frequencies the two signals are 90 degrees out of phase.

That is, for all positive $\omega$ that are overlapping, $h_1(\omega) = k_\omega h_2(\omega) e^{\pm i \tfrac{\pi}{2}}$, where $h_1$ and $h_2$ are the spectrums and $k_\omega$ is just a real scalar. The negative frequencies will be the conjugate of these.

If you mean correlation with signal shifts (cross-correlation - similar to convolution) then no, because if both signals have the same non-zero frequency then at some shift the sinusoids will not be 90 degrees out of phase, and therefore correlated.

More simply

Imagine a signal consisting of a sine wave and another consisting of a cosine wave of the same frequency. They will have overlapping spectrums (same frequency) but have zero correlation.

Add to the first signal a sinusoid of a certain frequency, and add to the second signal the same sinusoid but 90 degrees phase shifted.

Example

Barbara image with DC component removed, and Barbara image with DC removed and half spectrum multiplied by $i$ and the other half by $-i$ (90 degrees phase shift and conjugate for the negative frequencies)

Their spectrums overlap completely but they have zero correlation (at shift 0).


Why is a square wave aliased? I understand an ideal analog square wave contains sine-wave-components above Nyquist frequency (and towards infinity) and is thus subject to aliasing. So obviously we need an anti-aliasing-filter. Or do we? Why do we need one?

Intuition tells us that a square wave can hardly be aliased. It is so simple and there's no chance for a cycle or half cycle to be "skipped" unless the base frequency was above Nyquist... yet listening to a non low-passed-square-wave after sampling tells us otherwise. Many low pitches are heard, because of aliasing.

How is this possible?

Each successive sine wave seems to counteract the previous one, reducing the ringing at the nulls of the square-wave

Why does superposition hold on such a simple construct?

In addition to @hotpaw2 explanation, a graphic. There are two analog square waves (red and green), with different lengths. They are depicted with a fine sampling, denoted by crosses. Their actual sampling is denoted by circles. The  red one is shorter than the green one, as can be seen in the interval $]0.7\;0.8[$. Yet, the  sample points are the same. Thus, from the circled samples, you cannot tell one from the other. This  would happen for any square wave whose end falls inside the  blue segment. You can see here the loss of information.

With the  most simple schemes  (using cardinal sines), these different sampled squares would be reconstructed as a single analog signal (as their  discrete samples are the same.) Depending on the  precision on the sampling, several analog-digital-analog signals can be theoretically recovered,  possessing a certain amount of ripples, or overshoot, or Gibbs phenomena.

You can reduce the  ripples in the interior of the square wave (as seen going from blue to red to green), but the amplitude of the overshoot will remain constant  at the  vicinity of the  edges.

As commented by @robert bristow-johnson, there are a number of alternative techniques. The document Alias-Free Digital Synthesis of Classic Analog Waveforms (Stilson and Smith) provides some of them. I add a recent technique, based on Finite-Rate-Innovation: Feasibility of {FRI}-based square-wave reconstruction with quantization error and integrator noise, B. He et al., ICASSP 2015.

Finally, for your last question "Why does superposition hold on such a simple construct?" This can be regarded as a problem of approximating a non-continuous function on a basis on continuous ones: you can get pointwise convergence (all samples converge), but no uniform convergence (not all converge at the same rate).


Fast Cosine Transform via FFT I want to implement the Fast Cosine Transform. I read on wikipedia, that there is a fast version of the DCT which is similarly computed to the FFT. I tried to read the cited Makhoul* paper, for the FTPACK and FFTW implementations that are also used in Scipy, but I were not able to extract the actually algorithm. This is what I have so far:

FFT code:

def fft(x):
    if x.size ==1:
        return x
    N = x.size
    x0 = my_fft(x[0:N:2])
    x1 = my_fft(x[0+1:N:2])
    k = numpy.arange(N/2)
    e = numpy.exp(-2j*numpy.pi*k/N)
    l = x0 + x1 * e
    r = x0 - x1 * e
    return numpy.hstack([l,r])

DCT code:

def dct(x):
    k = 0
    N = x.size
    xk = numpy.zeros(N)
    for k in range(N):
        for n in range(N):
            xn = x[n]
            xk[k] += xn*numpy.cos(numpy.pi/N*(n+1/2.0)*k)
    return xk

FCT trial:

def my_fct(x):
    if x.size ==1:
        return x
    N = x.size
    x0 = my_fct(x[0:N:2]) # have to be set to zero?
    x1 = my_fct(x[0+1:N:2])
    k = numpy.arange(N/2)
    n = # ???
    c = numpy.cos(numpy.pi/N*(n+1/2.0)*k)
    l = x0 #???
    r = x0 #???
    return numpy.hstack([l,r])

*J. Makhoul, "A fast cosine transform in one and two dimensions," IEEE Trans. Acoust. Speech Sig. Proc. 28 (1), 27-34 (1980).

I've been reading about this and there are multiple ways to do it, using different size N.  My Matlab is rusty, so here they are in Python (N is length of input signal x, k is arange(N) = $[0, 1, 2, ..., N-1]$):

Type 2 DCT using 4N FFT and no shifts

Signal [a, b, c, d] becomes

[0, a, 0, b, 0,  c,  0,  d,  0,  d,  0,  c, 0, b, 0, a].

Then take the FFT to get the spectrum

[A, B, C, D, 0, -D, -C, -B, -A, -B, -C, -D, 0, D, C, B]

then throw away everything but the first [A, B, C, D], and you're done:

u = zeros(4 * N)
u[1:2*N:2] = x
u[2*N+1::2] = x[::-1]

U = fft(u)[:N]
return U.real

Type 2 DCT using 2N FFT mirrored (Makhoul)

[a, b, c, d] becomes [a, b, c, d, d, c, b, a].  Take the FFT of that to get
[A, B, C, D, 0, D*, C*, B*], then throw away everything but [A, B, C, D] and multiply it by $e^{-j \pi k \over 2 N}$ (half-sample shift) to get the DCT:

y = empty(2*N)
y[:N] = x
y[N:] = x[::-1]

Y = fft(y)[:N]

Y *= exp(-1j*pi*k/(2*N))
return Y.real

Type 2 DCT using 2N FFT padded (Makhoul)

[a, b, c, d] becomes [a, b, c, d, 0, 0, 0, 0].  Take the FFT of that to get
[A, B, C, D, E, D*, C*, B*], then throw away everything but [A, B, C, D] and multiply it by $2e^{-j \pi k \over 2 N}$ to get the DCT:

y = zeros(2*N)
y[:N] = x

Y = fft(y)[:N]

Y *= 2 * exp(-1j*pi*k/(2*N))
return Y.real

Type 2 DCT using N FFT (Makhoul)

Signal [a, b, c, d, e, f] becomes [a, c, e, f, d, b], then take the FFT to get [A, B, C, D, C*, B*].  Multiply by $2 e^{-j \pi k \over 2N}$ and then take the real part:

v = empty_like(x)
v[:(N-1)//2+1] = x[::2]

if N % 2: # odd length
    v[(N-1)//2+1:] = x[-2::-2]
else: # even length
    v[(N-1)//2+1:] = x[::-2]

V = fft(v)

V *= 2 * exp(-1j*pi*k/(2*N))
return V.real

On my machine, these are all roughly the same speed, since generating exp(-1j*pi*k/(2*N)) takes longer than the FFT.  :D

In [99]: timeit dct2_4nfft(a)
10 loops, best of 3: 23.6 ms per loop

In [100]: timeit dct2_2nfft_1(a)
10 loops, best of 3: 20.1 ms per loop

In [101]: timeit dct2_2nfft_2(a)
10 loops, best of 3: 20.8 ms per loop

In [102]: timeit dct2_nfft(a)
100 loops, best of 3: 16.4 ms per loop

In [103]: timeit scipy.fftpack.dct(a, 2)
100 loops, best of 3: 3 ms per loop


Can I decompose a sine signal into its component sine waves? If I have a sine signal that was composed of the sum of multiple other sine waves, is there any way that i can decompose the sine signal into its component sine waves without knowing any of the component sine waves or how many sine waves were summed to get the signal?
For example, if I create a signal be summing the following waves together: sin(x) + sin(3x) + sin(3x) + sin(4x), then can i take this signal, and without knowing the sine waves that were summed to create it or how many sine waves were summed to create it, can i decompose it into the sines that were used to create it (sin(x) + sin(3x) + sin(3x) + sin(4x) in this case)?

Curve fitting is another option after Fourier transform, you can fit a sum of sine function to your signal and the fitting coefficients show the amplitude and frequency of the signal.
Check this: https://stackoverflow.com/questions/38832474/curve-fitting-in-matlab-for-a-sinusoidal-function-with-more-than-8-terms


Guess stability based on filter coefficients is it possible to guess the stability of a filter with given coefficients a and b? I know that the right way would be to calculate the pole positions and check if they are in the unit circle but I think I read that there also is another method.

Cheers

In the general case there is no simple way to determine the stability of a discrete-time system (unless it is FIR, in which case it is always stable). In order to check stability you don't necessarily need to compute the poles of the transfer function, but you can use the Jury stability criterion, which does not explicitly compute the poles, because the latter is generally a numerically ill-conditioned problem.

For a second-order system the denominator coefficients have to be inside the stability triangle, which can be checked very easily.


This is how my professor is finding the frequency response of an LTI system when given the impulse response. Is this wrong? My DSP professor uses the freqz() function to find and plot the frequency response and I think it's wrong but I want to make sure before I point it out. Take this example:
Given the impulse response of an LTI, derive an expression for the magnitude and phase response and verify using MATLAB.
$$h(n)=2(-1)^{n-8\over 2}\sum_{k=0}^5 \delta(n-2k)$$
This is how I calculated the freq. response:
The impulse response is just the following values:
h = 2[1 0 -1 0 1 0 -1 0 1 0 -1]
The freq. response is the FT of the impulse response:
$$H(e^{j\omega}) = \sum_{n=-\infty}^\infty h(n)e^{-j\omega n}$$
h(n) exists only for $n=[0,10]$ so I change the limits of summation. Also, the series in the impulse response, converges to 6, so I replace it with 6:
$$\sum_{n=0}^{10} 2\times 6\times (-1)^{n-8\over 2}e^{-j\omega n}$$
Let $n'={n-8\over 2}\Rightarrow n=2n'+8$. Absorb the $(-1)^{n'}$ term. Replace these in the summation and update the limits in terms of n':
$$2\times 6 \sum_{n'=-4}^{1} (-1)^{n'} e^{-j\omega (2n'+8)}=12\sum_{n'=-4}^{1} e^{-j8\omega}(-e)^{-j2\omega n'}=12e^{-j8\omega}\sum_{n'=-4}^{1}(-e)^{-j2\omega n'}$$
And after finding the sum of the series, and some simplifications, the expression becomes:
$$H(e^{j\omega})=\frac{12(1-e^{-12\omega j})}{1+e^{-2\omega j}}$$
I derived the expression for phase response and got the following:
$$\angle H(e^{j\omega})=-5\omega +{\pi\over 2}+\left\{
\begin{array}{ll}
      0,   & k{\pi\over 6}\le\omega\le(k+1){\pi\over 6} \\
      \pi, & (k+1){\pi\over 6}<\omega<(k+2){\pi\over 6} \\
\end{array}
\right.$$ for $k=0,1,2,...$
I plug in some test value for $\omega$ to get the phase.
Then, to find and plot the frequency response, this is what my professor does:
freqz(h,1,w)
where h is the vector of the impulse response, h(n), and w is the vector of freq. (in rad). He says the '1' is added because it's an FIR filter.
Looking at the MATLAB documentation for freqz none of what he said makes sense. The '1' is supposed to be the arguments for the number of points for the transformation, and 'h' is supposed to be the digital filter represented either by the vector coefficients b and a or other methods.
When I run his code, I get the right shape and points where the frequency jumps up by $\pi$ but not the right values. For example, this is what I get with his plot:

Notice the phase at $\pi /6$.
This is what I get by passing the b and a coefficients from the transfer function $H(e^{j\omega })$:

Notice the phase at $\pi /6$. If I plug in $\omega = \pi /6$ in the expression for the phase response, I get that the phase is $\phi = -\pi /3 \Rightarrow -60^{\circ}$
I tried hours and hours to find out what my mistake was and got very frustrated, until I looked at the documentation. I want to bring this up to him because other students in the class are using his method which inherently gives wrong results, and future students will too, but just to make sure before I do that, am I right or am I wrong?
Thank you!

Your professor is right, and you're almost right too. The filter is clearly an FIR filter, but because its frequency response can be expressed as a geometric series, a recursive implementation is possible.
If you write the transfer function as a rational function you get
$$H(z)=2\frac{1-z^{-12}}{1+z^{-2}}\tag{1}$$
which is almost the same as you got, apart from a constant factor.
Note that there is a pole-zero cancellation in $(1)$, so the filter is FIR with impulse response
h = [2  0 -2  0  2  0 -2  0  2  0 -2];
So apart from rounding errors both

b = 2*[1,zeros(1,11),-1];
a = [1,0,1];
[H1,w] = freqz(b,a,1024);

and

[H2,w] = freqz(h,1,1024);

must give the same result, which turns out to be true. On my machine I get

max(abs(H1-H2))
ans =    7.9623e-14

Note that the phase is ambiguous, so the phase values are only identical up to an additive multiple of $2\pi$.


Calculating the peak power of a signal on python As far as I've researched, the energy and power of a given (discrete) signal are given by
$$E = \sum_n \left|x_n \right|^2$$
$$P = \lim_{N\rightarrow\infty}\frac{1}{2N+1}\sum_n \left|x_n \right|^2$$
Where N is the lenght of the given signal.
Im working with a signal embedded in some non-gaussian noise, and I want to calculate the ratio of the peak power of the signal and the power of the noise (see label of Fig 2 on https://arxiv.org/pdf/1701.00008.pdf).
Now, the power of the signal I already know how to compute, but I have not been able to find a precise definition of what peak power is, and I have not found a library on python to help me on that.
Is there a library on python that does this, or is there a concrete definition of peak power I can use to calculate it?
Thanks in advance!

Peak power is simply the highest square of any of the samples you encounter.
Average power is simply the average square.


Spectrogram of wav file I am calculating spectrogram of a audio file of 36 second using the following code snippet:

window = np.hamming(window_width)
length = signal.size
count = 0
for i in range(0, length-window_width, window_overlap):
    cur_frame = signal[i:i+window_width]
    window_frame = cur_frame*window
    spec = sp.fft(window_frame)
    if(i == 0):
        width = int(length/window_overlap)
        height = int(spec.size)
        specgram = np.empty([height, width])
        energy = np.empty(width)
    energy[count] = getEnergy(spec)
    for j in range(height):
        val = spec[j]
        specgram[j, count] = math.sqrt(math.pow(val.real, 2)+math.pow(val.imag, 2))
    count = count+1
print abs(specgram).max()

I get a ndarray of shape (4000L, 1008L)

but when i calculate spectrogram using Pylab:

Pxx, freqs, bins, im=pylab.specgram(sound_info, Fs=frame_rate)

I get ndarray of dimension: (129L, 125999L).

Actually I need to get the spectrogram which is the dimension of (frequency,nframes)

Thank you in advance.

How @Pichenettes said, we need more information to help properly

I'll take a guess:

If your Frame Rate are 44100, you will have 1587600 samples for 36 second of audio (36*44100), for the calculation example:

Pxx, freqs, bins, im=pylab.specgram(sound_info, Fs=frame_rate)

pylab.specgram seems use FFT Size = 256 and Overlap = 128 (50%) its give me the follow dimensional array:

floor(1587600 /  128) = 12403

(128, 12403)

You need set properly FFT Size and Overlap in your souce code.


How does `envelope` use the sliding window to determine RMS in MATLAB? I am trying to understand how does the sliding window works when using the function envelope in MATLAB to calculate RMS (root mean square). I'm using the following code to reproduce the behavior of envelope:

rng(1);
Array = randn(1000,1);
window = 100;
limit = length(Array)-window;

%Part 2
for i=1:limit
    RMS(i) = rms(Array(i:i+window-1));
end
plot(RMS)
hold

I expect the part 2 in the code above to generate the same results as:

upper = envelope(Array(1:limit), 100, 'rms');
plot(upper)

However the first option (blue line) is displaced about 50 samples compared to the second option (red line):

I would like to understand how envelope determines the first samples.

Since your window is not centered, you receive a shift of half-length of the chosen window.


State space representation in s-domain I was supposed to find state space representation and its matrices of this system:

and I have no idea, how to do this. We were told not to transfer the system to time domain, but I can only do state space representation from time domain schemas.

When I tried to solve this, I got matrices
$$A =  \left( \begin{array}{ccc@{\ }r}
    -a & k  \\
    -b & -p \\
    \end{array} \right)$$

$$B =  \left( \begin{array}{ccc@{\ }r}
    0  \\
    b \\
    \end{array} \right)$$

$$C =  \left( \begin{array}{ccc@{\ }r}
    1 & 0  \\
    \end{array} \right)$$

$$D =  \left( \begin{array}{ccc@{\ }r}
    0
    \end{array} \right)$$

Is that right? If not, how sould I solve it?

From visual inspection of the given block diagram, we obtain

$$\begin{bmatrix} s+a & -k\\ b & s+p\end{bmatrix} \begin{bmatrix} X_1 (s)\\ X_2 (s) \end{bmatrix} = \begin{bmatrix} 0\\ b\end{bmatrix} U (s)$$

which can be rewritten as follows

$$\left( s \begin{bmatrix} 1 & 0\\ 0 & 1\end{bmatrix} - \begin{bmatrix} -a & k\\ -b & -p\end{bmatrix} \right) \begin{bmatrix} X_1 (s)\\ X_2 (s) \end{bmatrix} = \begin{bmatrix} 0\\ b\end{bmatrix} U (s)$$

Hence,

$$\mathrm A := \begin{bmatrix} -a & k\\ -b & -p\end{bmatrix} \qquad\qquad \mathrm b := \begin{bmatrix} 0\\ b\end{bmatrix}$$

Since $Y (s) = X_1 (s)$,

$$\mathrm c^\top := \begin{bmatrix} 1 & 0\end{bmatrix} \qquad\qquad d := 0$$


how to detect and capture scene change in a recorded video I have a project to detect and capture scene change in a recorded video.Actually its a recorded video of a godown (of about 2 weeks) and I have to detect and capture pic whenever someone comes in.

I have tried vlc player motion detector. It detects changes very well but does not capture the detected scenes.
Is there any way that can help?

thanks and regards

You may use Low Rank Models for Video Background Subtraction  as a pre processing.
Once you leave out the background even simple methods will do to discern any motion.


How to apply "Substitution Method" in Difference Equations? How to get the value of $y_p(n)$ using substitution method?
What exactly do we need to substitute to get the value?

Substitution is not a method to obtain the particular solution $y_p(n)$, but the author simply claims that $y_p(n)=e^{-\beta n}u(n)$ is a particular solution of the recursion (1.70), and encourages you to check his claim by substituting $y_p(n)$ for $y(n)$ in Eq. (1.70) and see if it's correct or not:

$$e^{-\beta n}u(n)=e^{-\beta}e^{-\beta (n-1)}u(n-1)+\delta(n)=e^{-\beta n}u(n-1)+\delta(n)\tag{1}$$

Now check if (1) is satisfied. Distinguish 3 cases: $n<0$, $n=0$, and $n>0$. As a result you should be able to see that (1) (and, therefore, Eq. (1.70)) is satisfied for $y(n)=y_p(n)$, showing that the suggested sequence $y_p(n)$ indeed is a particular solution of (1.70).


How to filter a signal sampled at different frequencies I have a signal which was sampled using different frequencies. The 3 regions are:

0-10 ms: 10000 Hz
10-50 ms: 1000 Hz
50-200 ms: 100 Hz

The signal is very noisy, so I thought about smoothing it with a low pass filter. I worked with each of the 3 regions independently:

Apply window in time domain (Hann window)
FFT
Apply half Hann window, to remove high frequencies
IFFT

Each of my regions looks better now but I am having problems stitching them back together, due to the beginning and end of each stage deviating considerably from the original values.

What is the standard procedure of applying a low pass filter to a signal sampled at different frequencies during several stages?

Edit:

In case it might be relevant, this is the half Hann window used as a low-pass filter. The window length (in the example, 50) is different for every region. I am not interested in a specific cut-off frequency; the filter is used to remove noise and I am evaluating the output visually.

If I were you I would define my sampling period as the period present in the interval with the highest sampling frequency. Then I would try to interpolate the sections of the signal that were originally sampled at a lower sampling frequency to obtain a uniform signal of the same sampling frequency.

Then I would filter the whole signal with the same filter.

This procedure is often used when you have non-uniform sampling intervals but it might work for you.

As you say that the signal is very noisy maybe interpolating is not good but you can try to and see how it works.

If you can be 100% sure that decimating the intervals with higher sampling frequency is not going to introduce aliasing you can try to decimate all the intervals to the lowest sampling frequency and avoid interpolating noisy samples.

Maybe even an intermediate solution would be the best. Take every interval to the intermediate sampling frequency, decimating the one with higher sampling frequency (if you are sure there's no aliasing) and interpolating the one with lower sampling frequency.

Since your 3 sampling frequencies are integer multiples of each other you can easily decimate or interpolate from any one to get any other


Measure the length of curves via a smartphone's accelerometer Is it possible to measure the length of curves by moving a smartphone along them and processing the accelerometer's data?

For example, measure the circumference of a circle of radius $8 \,\rm{cm}$. I've read that accelerometers have noise but I've also read that their accelerometers are very accurate. What does it mean? Can we filter out that noise? What are the types of noise in this problem?

Accelerometers and gyroscopes provide corrupted measures of the acceleration $a$ and angular velocity $\omega$. So, they need to be preprocessed in order to use them.
Firstly, measure provided by the accelerometer ($a'$) is affected by Earth gravity $g$ (in the local frame reference), a bias $b_g$ and some random noise $\eta_a$. So if the true acceleration is called $a$, one has that measure is:
$$ a' = a + g + b_a + \eta_a$$
So in order to integrate accelerometer measures, one should know at least $b_a$ and $g$ in the accelerometer frame, and subtract them from $a'$. Moreover, in some datasheets you can find covariance of $\eta_a$. However, bias is time dependent term and it needs to be estimated each time we initialize the system, as well as the gravity in the local frame. In some cases, if the movement is not in gravity direction (2D movement) and gravity is aligned to some accelerometer axis (lets say $z$ axis), one can neglect $a_z$ and the knowledge of $g$ is not required (maybe it is your case).
Secondly, gyroscope is also corrupted by bias $b_g$ and noise $\eta_g$. So, if the true angular speed is $\omega$, one has that the measure is:
$$ \omega' = \omega +b_g + \eta_g$$
Again, you need to know at least the gyro bias. and subtract it from the measure.
Everything becomes easier if your device has no bias or it can be neglected w.r.t $a$ and $\omega$, and your movement is 2D (gravity is always oriented in the same way).
So once that you have your corrected measures, you can integrate them to obtain the trajectory:
$$ R(t+\Delta t) = R(t) \exp(\omega \Delta t)$$
$$ p(t+\Delta t) = p(t) +v(0) t + \frac{1}{2}R(t)a(t)t²$$
Where $R$ is the rotation matrix (see rotation matrix) and $\exp$ is the exponential map (see exponential map). $v(0)$ is the initial speed.
It can be concluded that having a known initial speed, a 2D movement, and a low corrupted gyro/accelerometer, it can be easily estimated the trajectory.
Edit: Calibration (i.e. getting $g$, $b_a$ and $b_g$) can be intermediately done if device remains static. In this case, accelerometer gives $a' = g + b_a +\eta_a$ and gyro gives $\omega' = b_g + \eta_g$