How do I warp a 1 dimensional signal in frequency domain? I want to implement VTLN (Vocal Tract Length Normalization).

It is equivalent to warping the signal in the frequency domain according to the following graph.

I can't wrap my head around how one implements such a transformation.

Suppose y = fft(x).

What should be y2 = f(y)?

What do I do after getting the warped frequency?

suppose x = [1, 2, 3, 4, 5]

y = f(X) = [ 10, 20, 30, 40, 50]

x_warped = [ 1, 1, 2, 4, 5].

What should y_warped[] be?

Do we delete some values and repeat some values in the y sequence?

is y_warped = [10, 10, 20, 40, 50]?

Image taken from

Pitz, Michael et al. “Vocal tract normalization equals linear transformation in cepstral space.” IEEE Transactions on Speech and Audio Processing 13 (2001): 930-944.

The mapping in this specific figure is basically a "lookup table" that maps $\omega$ to $\tilde{\omega}$ and the process in detail is outlined in section 5 of the paper that you link.

But on the point of applying that to a series, a function that effects such a transformation that is controlled by a parameter $\alpha$ is $f(x) = x^\alpha, x\in \left[0..1 \right]$ and $x \in \mathbb{R}$.

For example (Using Octave but easily extendable to other platforms):

function y=f(x,a) y=x.^a;end;
x = 0:0.001:1;
plot(x,f(x,1),x,f(x,2),x,f(x,1./2));
xlabel('x');ylabel('f(x)');
legend("x","x^2","x^{0.5}");
grid on;

Would give you:

Notice here that this is slightly different than the transform that you have with respect to the behaviour of the function with $\alpha$. I am using this for illustrative purposes here.

Notice also that because $x \in \left[0..1 \right]$, we can transpose that function to whatever limits we like. In your particular example it is $\pi$ or the number of bins of your Frequency Transform. Let's call that $N_{FT}$.

To apply the transform, you can either discretise the function or evaluate it on the fly. For example:

x_lookup = round(N_FT .* f(x/N_FT,alpha));

In this case, your input frequency transform bin would be x and that would have to now appear in bin x_lookup.

Notice here that x/N_FT to make sure that x stays within the range $0..1$ (assuming of course its lowest value is 0 here).

Again, section 5 contains all the details for specific f(x) functions you will need to effect this kind of mapping for a specific application.

Hope this helps.


FSK and IQ modulation I am having trouble understanding the FSK modulation process. There are many ways of demodulating the FSK signal.

What is the advantage of performing the FSK using IQ modulation?
How to demodulate the FSK signal?
Does the FSK signal look triangular wave or sinewave?

What is the advantage of performing the FSK using IQ modulation?

You only need one RF oscillator operating at a single frequency, instead of having 2 (or more in the case of M-ary FSK) oscillators operating at separate frequencies for each bit/symbol.
Since you only have one oscillator, you don't have to worry about discontinuities in the phase of the transmitted RF signal, like this:

Sharp discontinuities contain high frequency which widens the bandwidth of the signal. With I/Q modulation you'll have a much easier time producing a continuous-phase FSK signal (CPFSK).

  How to demodulate the FSK signal?

First lets assume we are trying to demodulate a simple FSK signal where the center frequency is $f_c$, a '1' is represented by the frequency $f_c+\Delta f$, and a '0' is represented by the frequency $f_c-\Delta f$.

Before IQ sampling, you tune your RF front end receiver to $f_c$. Then when you sample you have a baseband I/Q representation of the signal, centered around 0 Hz. The only thing to do is figure out if you have a positive frequency (+$\Delta f$) or negative frequency (-$\Delta f$). The easiest way to do that is to monitor the change in phase of the baseband IQ signal. A positive phase change corresponds to a positive frequency, and a negative phase change corresponds to a negative frequency.

So monitor the phase (the angle of $I + jQ$) and watch how it changes over one symbol period. Positive change = '1', negative change = '0'.

  Does the FSK signal look triangular wave or sinewave?

If you're talking about the RF signal, it looks like a sinewave of varying frequency.

If you're talking about baseband IQ signal, it also looks sinusoidal, but with a bunch of phase shifts, like this (blue = I, orange = Q):

In this particular example, a '1' corresponds to a phase change of $+\pi/2$, and a '0' corresponds to a phase change of $-\pi/2$.


Signal analysis: selection of oscillatory parts I'm trying to analyze the following signal:

My aim is to select the high-amlitude oscillation parts limited by the red stars.
At the beginning, i have tried to set a threshold and keep only the upper parts. however,a part that the threshold depend on the signal itself, this method didn't work well and it allows the selection of the non desired parts.
does anyone have another idea?

Thanks in advance.

By looking at that graph, it seems that simple math (or rather statistics) can be used to solve this problem. I'd divide (segment) your data into equal time frames and then compute the mean and possibly variance (or standard deviation) for each time frame.

The idea here is that the parts you're interested in will have a greater mean (and a greater variance) from the rest of the data. You'll still need a threshold, yes, but this approach is much more robust than just trying to compare the individual data points against the threshold.

For best results, you'll probably need to experiment a bit with different time frame lengths as well as different threshold values.


KLT tracker implementation I'm working on a object tracker based on KLT algorithm using matlab.
Can anybody guide me how I can find the features and calculate the distance in the sequence of the frames?
Thanks.

Try this example of tracking a single face using KLT.  Try this for tracking multiple faces.


RC circuit frequency response I am developing a project where I must analyze an incoming signal that was acquired from a microcontroller. The objective is to obtain the main frequency of the incoming signal.

At first, I’m analyzing an RC step response signal. So, basically, I want to obtain the frequency response, from the characteristic charging curve of an RC circuit.

My mentor told me that I could achieve this with the FFT. He told me that if I applied it twice at the incoming signal, I would obtain the frequency response curve. I think I did, but my only issue is the frequency axis. They’re completely messed up. Even when I apply the method to obtain the frequencies, it goes bad. (I know it's wrong because my cut-off frequency is around 11 Hz)

https://drive.google.com/file/d/1-ppjp6qE8J8O96Lj4uaJx36o2pt3KMSG/view?usp=sharing

So, what I’m asking is:

Why is it that applying the FFT twice to the signal give the frequency response?
Does anyone know how I can obtain the correct frequency axis?

Taking the FFT twice is similar to a method called "cepstrum".  It finds the spacing between harmonics.  You need to measure an oscillating signal to apply this.  I would go back to your mentor for more clarification.

Here is a solution for accurately measuring your charging rate in your diagram:  least square fitting to inverse exponential function

Ced

Followup:

I was very impressed by these articles.

From RC Charging Circuit in the section labeled: "RC Time Constant, Tau"

$$ V_c = V_s ( 1 - e^{-t/(RC)} ) $$

So when $ t = R C $ then $\frac{V_c}{ V_s} \approx .63 $

Since you already know $V_s$, the technique I referenced earlier is not necessary.  Looking at your graph $ R C \approx .07 $.

From Passive Low Pass Filter in the section labeled:  "Cut-off Frequency and Phase Shift"

$$ f_c = \frac{ 1 }{ 2 * \pi * R C } $$

Your cutoff frequency can then be found:  $ f_c \approx 2.27 $

From the section labeled:  "RC Low Pass Filter Circuit"

$$ X_c = \frac{1}{ 2 \pi f C } $$

$$ V_{out} = V_{in} \cdot \frac{X_c}{ \sqrt{ R^2 + X_c^2 } } $$

Divide the numerator and denominator by $X_c$:

$$ \frac{ V_{out} }{ V_{in} } = \frac{1}{ \sqrt{ \left( \frac{R}{X_c} \right)^2 + 1 } } $$

Substitute in $X_c$ and simplify:

$$ \frac{ V_{out} }{ V_{in} } = \frac{1}{ \sqrt{ \left( 2 \pi f R C \right)^2 + 1 } } $$

Now take the log (base 10):

$$ \log \left( \frac{ V_{out} }{ V_{in} } \right) = -\frac{1}{2} \log \left( \left( 2 \pi f R C \right)^2 + 1 \right) $$

From the section labeled: "Low Pass Filter Summary"

$$ Gain_{db} = 20 * \log \left( \frac{ V_{out} }{ V_{in} } \right) $$

Plug in the log of your voltage ratio to get:

$$ Gain_{db} = -10 \log \left( \left( 2 \pi f R C \right)^2 + 1 \right) $$

Your cutoff frequency is below the audible range of Hz, so you are going to attenuate all your frequency.  The higher the frequency, the greater the attenuation.  To find you fundamental frequency I recommend using a FFT and the frequency calculations I present in my blog articles.  You can find the link on my profile page.


What is the effect of sampling time, sampling frequency and SNR on latency? I believe that if we have two or more external sound cards, hence clock would vary for each device.

How will this affect the latency?
Will it increase or decrease?
How will the factors like

sampling rate,
sampling time,
SNR and
voltage reference impact latency?

...clock would vary for each device.

Unless you can afford high end studio devices that can synchronise to a common clock, this is taken for granted. However, saying that, you also need to balance your requirements. If you are using a relatively good "audio interface" (i.e. not a random off the shelf sound card) then you can expect a pretty accurate clock. Beyond this, the question would be: What is it that you do that requires audio synchronisation down to $\frac{1}{48kHz}$ seconds? I am taking 48kHz as a minimum for an audio recording here.

  How will this affect the latency?

The latency of each audio interface remains the same as before.

  Will it increase or decrease? How will the factors like sampling rate, sampling time, SNR and voltage reference impact latency?

For the latency to increase, you would have to introduce some sort of relationship between the audio interfaces. Instead, both are controlled by the same source and as far as each other is concerned, they operate in isolation. All they receive is a command ("start recording") and from there onwards, it is up to the computer/drivers/OS to be fast enough to serve the streams of data produced by each.

Sampling as far as each one of the interfaces is concerned is not affected. For the system overall, you might run into problems if you request the shortest buffer with the highest sampling rate but this then is up to the system the interfaces are attached to. For a given system you might have to try a few different settings to find the better one.

It is interesting that you mention SNR and voltage reference. The SNR of the signal path has nothing to do with the clock path, so, no influence there. Typically, the reference voltage (presumably of an ADC (?)) is inaccessible to you, so, again, no influence there.

If you are so worried that you might get some sort of phasing in the background because of the synchronisation, keep a "beeper" handy (like a clapper, but just for audio and even more accurate). After you hit "record", sound the beeper so that it is recorded through all channels (or at least through two groups that belong to the different interfaces you are trying to sync). Then, when the time comes to put everything together, just make sure that the beeper signals are aligned as best as possible in the beginning of the recording. From that point onwards, assuming good clocks for the individual interfaces, the recordings should remain "locked".

Hope this helps.


Time-shifting operation post the time-reversal operation when performing convolution I'm confused with the time-shifting operation post the time-reversal operation when performing convolution. Let's say we were to convolve $x(t)$ and $h(t)$, so I would have the term $x(k)$ and $h(t-k)$ to solve for. In $h(t-k)$, we first reverse $h(k)$ to $h(-k)$ and then shift it to $t$ units. Ideally, we should shift t units to the left but I find examples where they have moved $t$ units to the right instead.
I'm confused because, for any $x(t+n)$, we shift $n$ units to the left of $t$, so for $x(-t+n)$, why we are moving n units to the right of $-t$?

Just think of the value of $t$ where $x(0)$ appears. For $x(t+t_0)$ it is at $t=-t_0$, which corresponds to a left shift if $t_0>0$. However, for $x(-t+t_0)$ the value $x(0)$ occurs at $t=t_0$, which is to the right of its original position if $t_0>0$.
You can think of deriving $x(-t+t_0)$ from $x(t)$ in two different ways:

invert the time axis: $x(-t)$
replace $t$ by $t-t_0$, i.e., shift to the right if $t_0>0$: $x(-(t-t_0))=x(-t+t_0)$

Or

shift $x(t)$ to the left (for $t_0>0$): $x(t+t_0)$
invert the time axis: $x(-t+t_0)$: this flips the shifted function from the left to the right of the original function (again, if $t_0>0$).


Use of pulse shaping in digital modulation Sorry, I am a little bit confused here cause I thought when we say QPSK, the signal looks like the following.

I thought the bits to send/receive depends on the phase shift of the carrier, or is the above picture wrong or will never happen in real world?
So the discussion from here is to use pulses to send 1's and 0's? What does the signal look like in air actually? Pulse Shaped or like my attached picture? I am confused, so does the phase shifted carrier signal indicate a symbol or does the amplitude of a pulse shaped signal carries the bits?

QPSK in its raw form would appear just as the OP has shown, where the information is encoded into four symbols given by the four phase states alone. The problem is this has a very wide transmit spectrum given the pulse shapes are rectangular and the Fourier Transform of a rectangular pulse is a Sinc function (whose peaks roll-off in frequency only at the rate of 1/f). Wireless spectrum is an expensive commodity, so we go through great lengths in terrestrial wireless communications to limit spectral occupancy and in the process improve spectral efficiency. This is the only reason for pulse shaping, but it is significant. By more slowly transitioning between symbols rather than the abrupt transitions shown, we greatly limit the amount of spectrum needed (the spectrum will be given by the Fourier Transform of the pulse shape). Effective pulse shaping does this while controlling the inter-symbol interference that is added in the process.
The graphic below demonstrates this for a 16QAM waveform where the real portion of the time domain waveform is shown in the upper part of the graphic with and without pulse shaping. The lower part of the graphic shows comparative spectrum where we see the significant reduction in spectral occupancy due to the pulse shaping provided.

I see this graphic is the same one as what the OP linked, so to further explain: This is how the envelope of the signal would look like in the air (we don't see the actual carrier frequency in this graphic). This graphic was for QAM so does have information in the amplitude as well as phase, but for QPSK case the amplitude would carry no information. But even with the QAM case, the pulse shaping modifies the amplitude ONLY to reduce the spectrum, and modifies it in such a way to control the trajectory between our symbol samples of interest, while still passing through those samples defined by the modulation exactly at the proper symbol timing locations.
Below is an eye diagram showing the real and complex waveform for the QPSK case (The eye-diagram similar to the time domain waveform shown above, but repeating synchronous to symbol boundaries)

Show this shows how the over the air waveform would look (in actuality we typically do half the pulse shaping in the transmitter with a root-raised cosine filter while this plot shows the result after the second root-raised cosine filter in the receiver), but in general we are seeing the amplitude and phase of the time domain waveform for all samples, and the job of the receiver is to determine the sampling locations where the red dots are shown (the center of each symbol), and we then determine the magnitude and phase of those to demodulate (or demap) the QPSK data.
To see what is transmitted "over the air", the sketch below begins to demonstrate how the envelope of the carrier would be modulated by the pulse shaping, but doesn't provide the whole picture as the phase would also be transitioning from symbol to symbol so would not align with what would occur if we continued to transmit one symbol repeatedly. I added more details further below showing the exact waveforms expected. These are the real waveforms that are transmitted at the antenna rather than the baseband equivalent waveforms showing the complex envelope typically with two plots as I (real) and Q (imaginary). The sinusoidal carrier could be any frequency and the envelope would not change, including frequency = 0 which would then represent the eye diagram and I/Q plots that I have shown above (and must then be a complex waveform to represent).

A further note is the OP used symbol mapping that is not optimally gray coded (the 180 degree transition depicted in the diagram going from symbol 00 to 10 should be reserved for transitions where both bits change such that only one bit changes for all symbols that are close to each other).
Also note how this is similar to the operation of windowing prior to computing the FFT to reduce spectral leakage.

To avoid any confusion based on subsequent questions in the comments, here is the actual results on what the OP's QPSK waveform would look like at "RF" using properly gray-coded symbols using the following mapping of symbol to constellation:
0: -1-1j
1: 1-1j
2: -1+1j
3: 1+1j
We will start with the final waveform in the receiver as the desired QPSK constellation and work back to what is transmitted over the air by adding in the carrier offset and seeing the result:
First the I and Q waveform at baseband for a QPSK waveform pulse-shaped with a raised cosine filter with alpha = 0.3 would appear as the following in the receiver once all the carrier offsets are removed:

The constellation of the above is formed by plotting I vs Q on the complex plane. These trajectories represent the phase and amplitude of our RF carrier at every moment in time, and here I show the one instance that occurs at every integer symbol number in red, which is the only point over the duration of each symbol where the phase and magnitudes of the carrier will actually represent what we expect for QPSK (every other sample in between will be in transition at other phases and amplitudes!).

With an actual carrier frequency added, this same sequence would become the waveform shown in the plot below, where the amplitude of the envelope will match the amplitude of the constellation in the plot above:

Zooming in on the first 10 samples, this is what the OP would see for the sequence given (this is for the sequence 0,2,1,3,0,2,1,3,1,1,2,2.... specifically). The symbol sample points are aligned in this plot with the integer symbol locations on the horizontal axis.  This was done with a carrier frequency that is 10 times the symbol rate, but if we increased the carrier frequency further this plot showing the envelope would just get solid but appear identically the same otherwise.

And the same thing as we reduce the carrier to 4 times the symbol rate, showing that we will see the same basic envelope shape regardless of which carrier frequency is used (and as noted if the carrier gets too small compared to the symbol rate we need to use a complex signal to properly represent the waveform, exactly as we do in transceiver implementations):

Note that these are all done with a raised cosine filtered pulse-shape, but what is typically transmitted over the air is root-raised cosine filtered.  To be complete, below is the plot showing the same waveform with a root-raised cosine pulse shape:

With close observation of the above plot against the unmodulated carrier aligned with the first symbol we see how the modulated carrier is continuously changing in phase during the course of each symbol and, as also made clear by the constellation diagram of the modulated waveform, is only at the correct amplitude and phase according to the QPSK constellation only over very short durations, and in this specific case after being passed through a second RRC filter. (If the pulse-shape is altogether eliminated then the waveform will stay at the correct amplitude and phase for the entire symbol, and as the pulse shape roll-off factor is increased, this will approach that).

The above plots and descriptions used Raised Cosine (and Root-Raised Cosine) pulse shaping but an alternate harris-Moerder filter is described in this paper http://eon.sdsu.edu/~seshagir/SDR05.pdf. Using this as the pulse shaping filter results in an order of magnitude lower EVM for the same number of taps. This isn't commonly used even today (15 years later) to my knowledge since the root-raised cosine filter is baked into our specifications but certainly is very compelling to be aware of for consideration in future radios (I added a question to see if there is further experience with its use: Use of the harris-Moerder Nyquist Pulse Shaping Filter which may have additional information on this.)


How can I obtain the response signal for this question? In particular I am having trouble with 6b).

From what I understand, we can split a difference LTI equation into two sums, the sum of the previous responses, and the sum of the previous inputs.
Something like this:

Now to get a zero-state response, I have to zero out the first part of the formula. This I understand.

So to get the zero-state step response, I must use x[n] = u[n] and zero out all past outputs of the system.

The problem I'm having is deriving y[n] in the first place.

I have had a couple ideas, however am not confident in their correctness:

Since H(z) = Y(z)/X(z) I could compare the given H(z) to Y(z)(1 - 1/z) because that's the z-transform of 1/u[n]. Correct me if I'm wrong. Then I could make Y subject of the formula and figure out
the    inverse transform from there.
I thought about using convolution with the inverse transform of
H(z) and input signal u[n].

However I'm not convinced these are correct because nowhere am I "enforcing" the definition of zero-ing out the previous outputs.

Any guidelines would be appreciated!

HINT:

The (zero-state) step response is just the cumulative sum of the impulse response $h[n]$:

$$y_{ZS}[n]=u[n]\sum_{k=0}^nh[k]$$

This follows in a straightforward manner from the convolution of a unit step $u[n]$ with the impulse response:

$$y_{ZS}[n]=(h\star u)[n]$$


Why unit impulse function is used to find impulse response of an LTI system? Hello i am working in digital image restoration field, recently i have studied concept of convolution, i studied that to find the impulse response/point-spread function of an LTI system, an unit impulse function is passed from the system and its output is considered as impulse response, My question is that, Why an unit impulse function is used to find the impulse response of any LTI system, what is the mathematical and practical reason behind it.

It is important to realize that due to the LTI property the response to an impulse $\delta(n)$ describes the system completely, because any input signal can be written as a weighted sum of shifted impulses:

$$x(n)=\sum_kx(k)\delta(n-k)\tag{1}$$

Let $\mathcal{T}\{\cdot\}$ denote the system operator, so the response to a shifted impulse is

$$\mathcal{T}\{\delta(n-k)\}=h(n-k)\tag{2}$$

From (1) and (2), due the linearity and the time- (shift-)invariance of the system, the response to any input signal $x(n)$ can be written as a weighted sum of responses to shifted impulses:

$$y(n)=\mathcal{T}\{x(n)\}=\sum_kx(k)\mathcal{T}\{\delta(n-k)\}=
\sum_kx(k)h(n-k)$$

This last equality is the convolution sum describing the input-output relation of a discrete LTI system.


How to Estimate a Normalized PSD Using the Burg’s Method in MATLAB? I'll try to be most sintetic as I can.
I'm using an EEG system and I want perform an frequency analisys of the brain waves. I found this interesting article which describe how is possible to discern the "like" or "dislike" statement of mind looking at some specific band frequencies. To do this, they calculate the Power Density Spectrum (PSD) usign the Burg method and then normalize it. Following the relative extract from article:

  The power spectral density (PSD) of each epoch (16 epochs) for each
  subject (15 subjects) was estimated using the Burg method. The Burg
  method is an autoregressive (AR) pre-diction model based parametric
  spectral estimation method and yields a PSD estimate given by

  $\widehat P(f)=\frac{1}{f_s}\frac{\epsilon_c}{|1+\sum_{k=1}^{c}
\widehat a(k)e^{(2\pi jkf/f_s)}|}$

  where fs is the sampling frequency, c is the order of the model, f is
  the frequency, and a(k) are the AR model parameters. These parameters
  are computed by the minimization of the back-ward and forward
  prediction errors while satisfying a specific recursion scheme called
  the Levinson–Durbin recursion. In the Burg method the signal whose PSD
  to be estimated is assumedto be the output of an all-pole linear
  system driven by whitenoise. In our study, the model order was taken
  as 15 after literature. This method estimates the model
  parameters directly without a need for the autocorrelation function
  calculation. The main advantages of this approach are
  resolving closely spaced sinusoids in signals with low noise levels,
  and successfully estimating power spectral densities of short
  data recordings. In addition, this technique guarantees a
  stable autoregressive model and is computationally efficient. For each
  epoch, we computed/estimated PSD values at each integer frequency. To
  be more specific, we have computed thespectrum at f = 4, 5, 6, . . .,
  39, 40 Hz, total of 37 frequency points. We then normalized the power
  values using the sum of all power values in that segment to minimize
  inter-subject and intra-subject variability. Therefore, we obtained 37
  normal-ized power values (NPVs) for each epoch (totally 16 epochs),
  [...]

Now I want recreate this analisys in MATLAB but I haven't deep knowledge of this method. I found only an specific function but I don't know how to use it to have the normalized values (NPV's) described above. I tried it but the I don't understand for example why the result band frequency is 0-129Hz.
The lenght of my epochs is 1000 sample whit a sample frequency of 250Hz.

Maybe the answer to my question could be long and complicated but can anyone explain me how to use the MATLAB funcion correctly? Why in the article choose the Burg's method to obtain the PSD? And what means

  We then normalized the power values using the sum of all power values
  in that segment

?

The link to the paper is not working.

If you could show your code maybe we could help to correct it.

With Fs=250Hz you will only be able to see frequencies up to 125Hz (Fourier transform theory).

For using the pburg method you only need to provide the samples and the order of the model (usually the model order is the double of the number of spectral peaks).


Why is my fundamental frequency at the 15th harmonic order I am trying to implement a code in MATLAB to show harmonic content and order in Matlab. my fundamental is supposed to be at harmonic order 1, however the ouput plot I am getting is wrong and I am having difficulies proceeding.
I am a first time user of this forum and i have posted the same question in the matlab site
Here :
I would really apreciate any help. Thanks

[D,S,R] = xlsread('test data.xls');
v = D(:,2);
Signal = D(:,2);
Ts = 0.00005;                                               % Sampling Interval (seconds)
Fs = 1/Ts;                                                  % Sampling Frequency (Hz)
Fn = Fs/2;
%PERFOM FFT
N = length(Signal);
meanSignal = mean(Signal);                                  % ‘Signal’ Mean
FTSignal = fft(Signal-meanSignal)/N;
Fv = linspace(0, 1, fix(numel(FTSignal)/2)+1)*Fn;           % Frequency Vector
Iv = 1:numel(Fv);                                           % Index Vector
%Plotting Harmonics (THD)
har_mag = abs(FTSignal(Iv))*2;
figure;
bar(har_mag(1:20),0.4)
xlabel('Harmonic Order','fontsize',10);
ylabel('Voltage (V)','fontsize',10)
hold off

False reading of the output I am getting is this:

the output should look similar to  figure below:

Yes you are correct that you are indeed seeing the equivalent of what would be the 15th harmonic of the fundamental frequency of the equivalent Fourier Series Expansion of the resulting equivalent continuous time waveform that was truncated in time to 40 samples (2 ms). According to the Fourier Series Expansion, $T$ would be 2 ms and the fundamental frequency for any continuous time domain trajectory over duration $T$ would be $1/T$. That fundamental frequency need not be the largest frequency nor need it be non-zero. The sum of that frequency with it's initial phase and amplitude with all the higher integer harmonics will converge to the original waveform.

That is a rather convoluted but mathematically accurate view. You are also seeing the expected result of the Discrete Fourier Transform for what must be a dominantly sinusoidal waveform that has a frequency of 15/40 of the sample rate (at or near 7.5 KHz). You are not sampling high enough to uniquely see higher order harmonics as they would be folded into the first Nyquist Zone from DC to $f_s/2$. (For example the 2nd harmonic at or near 15KHz would fold to 5KHz given your 20KHz sampling rate, or bin #10 in the plot above, but the signal may also have phase or amplitude variation causing such a sideband that is 2.5 KHz away- once sampled you can't distinguish these two effects).  You could certainly rotate (frequency translate) your signal to move to a lower bin as in your desired plot, but that will not resolve the ambiguity of the harmonics unless you do that frequency translation prior to sampling (as well as low pass filter out any signals that may still exist above $f_s/2$.)


How to shift the frequency spectrum? Lets say we have a spectrum ranging from -X MHz to +X MHz. I would need to correct the frequency error in the spectrum by shifting the zero component to the middle (0 Hz).

If the output (the frequency spectrum) is calculated via FFT, as far as I know I can move the spectrum by adjusting the 'twiddle factors' (or coefficients, for complex data sine and cosine waves).

In the case of a size 1024 FFT (bin indexes from 0 to 1023), 0 Hz component should exist in bin number 511. However, due to possible frequency error the 0 hz component may actually be in bin 510 for example.

I cannot seem to find much information on this. Any help appreciated.

EDIT: Mistake in the question.

If the frequency shift that you want is a multiple of the bin spacing, as in your example, then you can easily effect the shift that you want by just rotating the FFT outputs by the number of bins that you need. In the more common case that the frequency offset is not an integer multiple of the bin spacing, then you can multiply the signal by a complex exponential function before doing the FFT.

So, if you determine that the center frequency component that you speak of is actually located at frequency $f_{offset}$ Hz in your data, and the data is sampled at rate $f_s$ Hz, then to shift the spectrum such that the component of interest is at zero frequency in the FFT output, you would do:

$$
x_{\text{shifted}}[n] = x[n] e^{-\frac{j 2 \pi f_{\text{offset}} n}{f_s}}, \ \ n = 0, 1, \ldots , N-1
$$

$$
Y_{\text{shifted}}[k] = FFT\left[x_{\text{shifted}}[n]\right]
$$


Time delay estimation between two complex signals(IQ) I need to time align two complex signals. The two signals obtained were at different RF frequency before digitisation. So,I have two complex signals(IQ) and I need to do delay estimation and apply correction on one of the signal.

Cross-correlating the absolute values of two signals and finding the peak gives one answer while cross-correlating the complex values and then finding the peak from the absolute value of the obtained (complex)correlation gives another peak?

Which method shall I choose? what is the best way forward? Any insight into it?

Cross-correlating the absolute values of two signals will often not yield good results.  Consider the case in which you are trying to align two continuous (i.e., non-pulsed) signals.  Their absolute values are roughly constant with time.  Cross-correlating them will yield no peak.  If the amplitude of the signals changes with time, you may get a peak, but the cross-correlation will be able to exploit only the amplitude of the signals.

On the other hand, cross-correlating the complex values will yield a better result because it exploits both the amplitude and phase of the signals.  Note that you will generally need to look at the absolute value of the correlation, not just the real component, since the phases of the two signals are usually not the same.  Some complication arises if the two time-delayed copies of the signal have different frequency shifts, in which case long correlations may not yield a peak.  That frequency difference can sometimes yield useful information, however.  This is a significant topic in radar.


Sharpen Defocused Image (Deconvolution / Image Restoration) Using OCR, I want to extract text from product packages using Google Glass. However, because of the fixed focus of the camera the package pictures are blurred. Is there a way to sharpen the image? Currently I use unsharp masking to enlarge the gradients of the edges, which gives me OK results.

Is there a better way to do this? I thought about taking a picture of a point-spread-function and using this to deconvolve, but I doubt this gives me a good approximation of the distortion kernel.

Here is an example:

You're after an algorithm in the family of "DeConvolution".
Specifically in your case, is called Blind Deconvolution.

Yet if you have some assumption the Blur you can use Wiener Filter or Lucy Richardson.
Both of them are actually the MMSE Estimator just with different assumption of the noise.
Both of the methods are actually "Inverse Filter" on an Low Pass Filter, which means they are High Pass Filter, just like the Unsharp Filter you applied.
The difference is those methods are optimal per given "Blur Model".

Yet, if you go after the Blind Deconvolution, today the best estimators are based on some prior of the image.
Something like the distribution of the Gradient of the image and stuff like that.

Yet if you are only after the text in the image.
Something simple like Wiener Filter + Edge Detection should do the work to give you most data.


Questions about reproducing results of a paper on Communication Systems and Wireless Communication I'm trying to reproduce the results of a paper on Device-to-Device Communication, channel assignment and mode selection (G. Yu, L. Xu, D. Feng, R. Yin, G. Y. Li and Y. Jiang, "Joint Mode Selection and Resource Allocation for Device-to-Device Communications," Nov. 2014.).

I have a couple of questions on the noise spectral density, how to generate the channel taps, and finally the throughput of the entire system.

It is mentioned the following in the paper:

  "All links are assumed to experience independent block fading.
  Hence, the instantaneous channel gain of the interference
  link between CU m and the receiver of D2D pair k can be
  expressed as
  $$h_{k,m} = Gβ_{k,m}d^{−\alpha}_{k, m},$$
  where G is the path loss constant, α is the path loss exponent,
  $β_{k,m}$ denotes the channel fading component, and $d_{k,m}$ is the
  distance between CU m and the receiver of D2D pair k.
  Similarly, we can express the channel gain between CU m and
  the eNB as $h^C_{m,B}$, the channel gain between D2D pair k as $h^D_k$,
  and the channel gain between the transmitter of D2D pair k and
  the eNB as $h^D_{k,B}$."

Towards the last section, we have the following relevant simulation parameters:

Noise Spectral Density: -174 dBm / Hz.
Path Loss Model (for Cellular Links): $128.1 + 37.6 \times \log_{10}(d[km])$.
Log Normal Shadowing: 10 dB.
Uplink Bandwidth: 3 MHz.

Now my questions are related to the channel model:

How do I convert from the Path Loss Model combined with the Shadowing Effect of 10 dB to the first mentioned equation ($h_{k,m} = Gβ_{k,m}d^{−\alpha}_{k, m}$) ?
When I convert the noise spectral density from dBm/Hz to dBm, should
I divide the value by the number of Uplink channels, which is 20 ?
And if so, should I do that in the dB or Watts units ?

Let's tackle this piece by piece:

Path loss

Your path loss model gives the path loss in dB. Convert to linear scale. You can directly read $\alpha$ and $G$ from the result.

Haven't read that paper (paywalled), but it looks to me like your shadowing should be part of the fading  $\beta$. Make sure your channel model is compatible with the one used in the paper at all! There's no guarantee that different sources use channel models that are.

Power Spectral Density to Power

Well, assuming white noise (check that assumption is correct for this paper!), the noise power in a channel is of course bandwidth of that channel times noise power spectral density.

  should I do that in the dB or Watts units ?

dB is not a unit!

Generally, the actual units don't matter, as long as you do all the calculations correctly, so it really doesn't matter whether you're working in W or (MeV per minute).

Your first question, however, points to you not being 100% familiar with what dB and dBm are. I strongly recommend reading the Wikipedia article on decibel if that is the case!


Proving for windowing function Prove that $1+2\cos(\omega)$ is the same as $\frac{\sin(\frac{3\omega}{2})}{\sin(\frac{\omega}{2})}$

Knowing that
$$p[n]=\sum_{k=-1}^1 {\delta[n-k]}=\delta[n-1]+\delta[n]+\delta[n+1]$$

Doing the FT yields
$$\begin{align} P(e^{j\omega}) &=\sum_{k=-\infty}^\infty p[n]e^{-j\omega n} \\
&=e^{-j\omega}+1+e^{j\omega} \\&=1+2\cos(\omega) \end{align}$$

But, how does this fit into the general formula of $\frac{\sin(\omega(M+\frac{1}{2}))}{sin(\frac{\omega}{2})}$ where, in this case my $M=1$?

Thank you.

I have come across my ways.

$$\begin{align} P(e^{j\omega})&=1+2\cos(\omega) \\ &=\sum_{k=-\infty}^\infty \sum_{k=-1}^1 \delta[n-k]e^{-j\omega n} \end{align}$$

Knowing that non-zero only when $n=k$,

$$P(e^{j\omega})=\sum_{k=-1}^1 e^{-j\omega n}$$

Let $M=n+1$,

$$\require{cancel} \begin{align} P(e^{j\omega}) &=\sum_{k=0}^2 e^{-j\omega (M-1)} \\ &=e^{-j\omega (-1)} \sum_{k=0}^2 e^{-j\omega M} \\ &= e^{j\omega} \left(\frac {1-e^{-j\omega 3}} {1-e^{-j\omega}}\right )  \\&=\cancel{e^{j\omega}} \left( \frac {e^{j\omega 3/2}} {e^{-j \omega 3/2}} \right) \left( \frac {e^{j\omega 3/2} - e^{-j\omega 3/2}} {\cancel{e^{j\omega}} e^{j\omega 1/2} - e^{-j\omega 1/2} } \right) \\&= \frac{\sin(\frac{3\omega}{2})}{\sin(\frac{\omega}{2})} \end{align}  $$


how can I calculate pitch deviation in a large collection of audio files? I am making a software model of an existing digital synthesiser and I have made audio recordings of single notes with particular parameter settings. I made a program in C++ which sends a MIDI sysex message to change the parameter in question, then triggers a MIDI noteon and noteoff message. The sound is recorded to a wav file. So I have several collections of 100 wav files for each parameter (the parameters in question have 100 steps).

One of the (easier) parameters I am looking at is LFO depth. The LFO is applied to the oscillator pitch. So I have 100 wav files with a square wave LFO modulating the oscillator pitch between a base pitch and another pitch that goes between base pitch (zero modulation) and base pitch+1 octave. I would like to write a script to work out precisely how many semitones/cents modulation occurs in each wav file. So I would like to end up with 100 numbers that I can use to map to modulation depth in my software model.

Can anyone recommend a tool/and or libraries for doing this task?

Use a very slow LFO to get very long segments of uniform pitch (this will allow you to use long analysis windows and counter 'border' effects), and use the simplest possible waveform for the oscillator to prevent octave errors in pitch detection. Aubio provides an implementation of the YIN algorithm which is fairly accurate.

I personally doubt there's anything fancy going on here - it is very likely that the relationship between the LFO modulation depth and the pitch deviation in cents is linear.


Amount of exactly approximated samples with Padé approximation I've already had some good discussion with Fat32 in this question yesterday.

Today I'm confused again. The author of Statistical digital signal processing and modeling M.h. Hayes and Fat32 both stated that Padé can approximate p+q+1 samples. Then I found this exercise

and the solution

Do I count wrong or are there 5 samples used for the approximation? q+p+1 should be $2+3+1=6$ and therefor 6 perfectly approximated samples. Okay I've calculated it. Due to the periodicity the sixth sample is correct too. But what if the vector would be e.g. 10 at $x[5]$ (sixth sample)?

To solve for the a's I'd then have 3 equations for the vector $[1, a_1, a_2]^T$ which is impossible to solve. Where do I make the mistake?

The problem and its solution is consistent with the Monson Hayes' book that describes Padé approximation.

Your mistake is that you've taken $p=2$ but $q=3$ which is wrong. Correct is $q=2$ and there are 3 unknowns in $b[k]$ , that the problem defines a second order system , i.e, $a[k] = [1, a_1, a_2]$ and $b[k] = [b_0, b_1, b_2]$ , according to the book.

Since there are $p+q+1 = 2+2+1 = 5$ (that +1 is for $b_0$) unknown coefficients in the system transfer function $H(z)$, then from its impulse response $h[n],$ only the first 5 samples will match to the data $x[n]$ for $n = 0,1,2,3,4$.

And from the given solution as you can see there are 5 equations for 5 unknowns, of last 2 is used to solve for the unknowns $a_1,a_2$ and the first 3 used to find $b_0,b_1,b_2$. The data set used is the $x = [x_0, x_1,x_2,x_3,x_4] = [2,1,0,-1,0]$.

Note that for the first two equations $x(-1)$ and $x(-2)$ are implictly referenced, which are assumed to be zero for the causal modeling.


Discretization of Sinusoidal Signals How do I should sample two continuous sinusoids of $f_1$=1/4 Hz and $f_2$=1/2 Hz at a 1 Hz rate?

Can i recover both if i bandlimit after sampling at a 1/2 Hz limit?

I have:

The function $h_1(t)=\cos(2\pi\frac14 t)$ sampled by $|||(t)$ has fourier transform $H_1(f)=4||(4f)*|||(f)=2|||(2(f-1/4))$.
The function $h_2(t)=\cos(2\pi\frac12 t)$ sampled by $|||(t)$ has fourier transform $H_2(f)=2||(2f)*|||(f)=|||(f-1/2)$,

When I bandlimit $h_1$ with $\Pi(2f)$ I obtain $4||(4f)$.
If I try to bandlimit $h_2$ by the same window, exploiting the Heaviside half maximum convention $\Pi(\pm1/2)=1/2$, I obtain $2||(2f)$.

Both sinuoids are perfectly recovered, though the Heaviside limit in the second case is not applicable in practice for dozens of reasons.

In here, I am using the pulse, sampling and sinusoidal symbols. $H$ is the Heaviside unit step function:
$$
\Pi(f)=H(f+1/2)-H(f-1/2),\\
|||(t)=\sum_{t_n \in \mathbb{Z}}\delta(t-t_n)\\
||(t)=\frac12\delta(t+1)+\frac12\delta(t-1)
$$


Calculating spatial width from power spectrum HWHM Say I have a 1D (spatial) signal (resolution = $1000$) which is zero everywhere except from $x = 250$ to $750$, where it equals one.

I ultimately want to calculate the spatial width of this signal using FFTs. Of course we know the width here to be $500$; in actuality, I am dealing with a signal that evolves with time and wish to calculate the average "pulse" width over all the time frames, so I do not know the widths. I have opted to use FFTs in this pursuit, so I must conduct a "sanity check" to make sure the method works. This method was suggested to me by a colleague whose intuition is many leagues farther than my own, so if someone could explain the intuition to me, I would appreciate it a lot.

Step 1: Subtract the DC background (subtract the mean from every point of the signal).
Step 2: Take the FFT of the signal, then the power (the Fourier transform times the complex conjugate of it). Normalize the power spectrum.
Step 3: Calculate the half-width at half-maximum (HWHM); here half-width is the half-width of the peak in k-space, of course.
Step 4: Convert this k-space HWHM back to real-space: real-space width = 1 / (HWHM / resolution).

When I do these steps for the signal above, I calculate a real-space width of $1189427$, laughably off from $500$. Where does the method go wrong?

Say I have a 1D (spatial) signal (resolution = $1000$) which is zero everywhere except from $x=250$ to $750$, where it equals one.

This is not "resolution". Resolution is 300 Dots Per Inch. In which case, we could say that the total physical length of your pulse is $\frac{500}{300} \approx 1.666$ inches (or any other unit of length).

  I ultimately want to calculate the spatial width of this signal using FFTs.

Why?

  I am dealing with a signal that evolves with time and wish to calculate the average "pulse" width over all the time frames, so I do not know the widths.

If there will be multiple pulses of different widths on the same signal, then by opting to detect them with the FFT you are setting yourself a very big challenge because the FFT would return to you information about the signal as a whole. So you could, for example, derive an average rate of pulses (even using the algorithm that is presented here) but not the widths of individual pulses.

If it is somehow guaranteed that within a window of 1000 samples, there will be a pulse whose length is guaranteed to be staying well below 1000 samples and all we have to do now is detect where the pulse is and how long it is, then opting for the FFT is an overkill.

The usual way to detect pulse widths is via the simple use of a threshold and a counter. Once the signal's amplitude goes above the threshold, the counter starts counting and it stops once the signal's amplitude goes below the threshold. If you are going to operate in a noisy environment, then there are a number of improvements to that such as adding hysterisis to the threshold, so that it doesn't respond to very short "bounces" of the waveform and adaptive thresholding where the threshold limit would be derived from the given window of observation (here, from the 1000 samples).

If you absolutely have to work in the frequency domain, it might be better to look into the discrete wavelet transform (DWT), whose output is a time/scale(frequency) representation. But the actual detection of the pulse width is likely to be happening (again) using some form of threhsolding on the output of the DWT. (So, again, huge overkill).


Why does my digital upconversion result in two tones? I am using Xilinx's RFSoC in I/Q mode. The important point is that in the digital space I have real and complex samples. I take the incoming data and multiply it by a complex sinusoid.

In my test setup I am sampling a 125 MHz sine wave and am mixing it with a 5 MHz sinusoid. I take the mixed signal and pass it out a DAC to view it on a spectrum analyzer. What I would like to see is the 125 MHz digitally upconverted to 130 MHz. What I actually see is two tones: 120 MHz and 130 MHz.

Is there a way, without just filtering the 120 MHz sine wave, to upconvert the 125 MHz tone to 130 MHz by mixing?

You can use the property of Fourier transforms which says that: $x(t)e^{j2\pi f_0t}\leftrightarrow X(f-f_0)$. That is, multiplying the time domain signal with a complex exponential will give you a shift in the frequency domain.

By multiplying by another sine, you get: $x(t)\text{sin}(2\pi f_0 t) \leftrightarrow X(f)*\text{FT}\big(\text{sin}(2\pi f_0 t)\big)=\frac{X(f-f_0)-X(f+f_0)}{2j}$, and so you see an extra peak at $120$ MHz.


Is autospectrum the same as power spectrum density? In my slides about signal processing, there is one that mentions the same thing as the beginning of this and this answer, namely that the Fourier transform of a signal, squared is the power spectrum density of the signal.

In this talk, it is mentioned that coherence is calculated by dividing the squared cross-spectrum by the product of the two auto-spectra.

However, the formula in my slides divides the squared cross-spectrum by the product of the formula of our previously seen PSD and another PSD.

So, is autospectrum the same as PSD? I can find a lot of information about PSD, but not on autospectrum.

This may be just a matter of semantics. Indeed it is true that "The Fourier Transform of the auto correlation x[n] is identical to the magnitude squared of the Fourier Transform of x[n]". This is just a mathematical identity.

You could call that the power spectrum density (PSD) but in most practical applications any actual PSD computation will involve some sort of framing and windowing. The choice of framing and windowing will impact the result (and it's a complicated trade off), so there isn't really one clear and unambiguous definition of PSD.

You can still use the mathematical identity but it needs to be properly adjust with respect to framing, windowing and circular vs. linear convolution.


How to reconstruct RF signal using IQ data I have a complex IQ signal: $s =  X  + Yi$. I know my original $f_0$ was $5\textrm{ MHz}$.

I can calculate the magnitude/amplitude of $s$: $A = \sqrt{X^2 + Y^2}$
and phase, $\phi = \arctan\left(\frac YX\right) $

But how can I reconstruct the signal using this frequency, amplitude and phase data?

I need the output as a time varying signal which incorporates the amplitude and phase information from the IQ data. I have found this:

  The reconstruction of RF-data from IQ-data is straightforward. It is a
  reversal of the complex demodulation  in  the  previous  section.  The
  decimation  is  reversed  by  interpolation.  The  low- pass  filter
  cannot  be  reversed,  but  should  be  chosen  without  loss  of
  information  in  the  first place.  The  down  mixing  is  reversed
  by  up  mixing.  At  last,  the  RF-signal  is  found  by  taking the
  real-value of the complex up-mixed signal.

But I could really do with a worked example.

EDIT:

This MATLAB code appears to do what I want:

IQ = resample(IQ,downsampl,1);

[Z,X] = size(IQ);
time = [0:Z-1]'*1/F * ones(1,X);

signal = IQ .* exp(+2*pi*i*time*carrier);

signal = real(signal);

where resample resamples IQ at downsampl/1 times.

It is not clear to me what I should use for downsampl and what F is?

The RF signal $r(t)$ is obtained from the complex baseband (IQ) signal $s(t)=x(t)+jy(t)$ in the following way:

$$r(t)=\text{Re}\{s(t)e^{j\omega_0t}\}\tag{1}$$

where $\omega_0$ is the carrier frequency (in rad/s). There are two other equivalent representations of $(1)$, where I use $s(t)=x(t)+jy(t)=a(t)e^{j\phi(t)}$:

$$\begin{align}r(t)&=x(t)\cos(\omega_0t)-y(t)\sin(\omega_0t)\tag{2}\\
r(t)&=a(t)\cos(\omega_0t+\phi(t))\tag{3}\end{align}$$

Eq. $(2)$ shows that the I- and Q-components are both modulated by orthogonal carriers, whereas $(3)$ shows that $r(t)$ generally exhibits amplitude modulation as well as phase modulation.

The quote in your question suggests that in the respective text the signals are all represented as discrete-time signals, hence the decimation and interpolation stages.


How can we show that $|H(e^{j\omega})|=|H(e^{-j\omega})|$? Let $H(z)$ be the rational system function of an LTI system.
How can we show that $|H(e^{j\omega})|=|H(e^{-j\omega})|$?

Your statement only holds true for a real sequence $h[n]$.
The frequency response of $h[n]$ equals to
$$H(e^{j\omega}) = \sum_{n=-\infty}^{\infty}h[n] e^{-j\omega n}$$
and if $h(n)$ is real, the conjugate symmetry condition holds
$$H(e^{-j\omega}) = \sum_{n=-\infty}^{\infty}h[n] e^{j\omega n}=\big(H(e^{j\omega})\big)^*$$
So we can derive that
$$\big|H(e^{j\omega})\big| = \big|H(e^{-j\omega})\big|$$
and
$$\arg\{H(e^{j\omega})\} = -\arg\{H(e^{-j\omega})\}$$
where $\arg\{\cdot\}$ represents phase response.
We can also derive that the real part of the frequency response is even symmetric and the imaginary part is odd symmetric.
$$\Re\{H(e^{j\omega})\} = \big|H(e^{j\omega})\big| \cos\big(\arg\{H(e^{j\omega})\}\big) = \Re\{H(e^{-j\omega})\}$$
$$\Im\{H(e^{j\omega})\} = \big|H(e^{j\omega})\big| \sin\big(\arg\{H(e^{j\omega})\}\big) = -\Im\{H(e^{-j\omega})\}$$


How can the FFT be used for estimating linear chirp parameters? Consider a linear chirp signal, i.e. $f(t) = f_0 + kt$, where $f_0$ is the starting frequency, and $k$ is the chirp rate. After applying an FFT, how would I find the starting frequency and chirp rate?

Why use an FFT for this? I think there are better simpler ways to do this. For example, you can simply measure the distance between zero crossings to get a first order estimate for frequency vs time and then refine with reconstructing the sweep and looking at the difference of the original and the reconstructed signal.


If I can make a call to another city, why does wifi only reach 50 meters? What makes it possible for me to communicate with someone in another city when my phone can only reach the modem / wifi networks that are within 50 meters?

The signal may be sent between multiple destinations before reaching the final destination / end point, but I assume the first destination is still more than 50 meters away.

This question, having a physical and network rather than mathematical nature, is close to being off-topic, but for your convenience let me put an answer for it.

Making a phone call is traditionally a communications network operation. But recently IP telephony over computer networks (internet) also emerged and furthermore GSM (mobile cellular) systems are so much intermingled with internet that the traditional medium seems to have blurred boundaries.

Making a traditional phone call using your GSM cellular phone, involves accessing an operator base station (BS) closest to your location. Base Stations, work at 900/1800 MHz, ranges and can be reached at as low as about -110 dB signal levels. Therefore, your phone can reach a BS about a 5/10 kilometers away in free space and less than a km in city areas. Once it's reached, remaning connection to other cities etc., is powered and handled by the BS network itself. Your phone is only required to reach the closest BS.

Wi-Fi is term used to denote wireless internet access point, that enables internet capable devices to wirelessly connect to internet. It may be at your home, office or airport, operating at 2.4 (or recently 5) GHz. It's signal output power level is low, therefore it can be reached within a hundred or less meters for fast / reliable connections. Once you access a wifi, the remaining communication is powered by Wi-Fi adapter / router device and your device is not involved.

As another option, you can also make an IP phone call using WhatsApp, Skype etc with your smart phone-computer combo. Then it's either accesing the closest (and allowed) Wi-Fi access-point at home/office or it may also be doing it through the GSM BS internet using its mobile data feature...


Bit-reversal equivalence on IFFT (radix-2 Cooley-Tukey) From everything that I read on the internet, I understood that when performing FFT or IFFT, one should perform bit-reversal to obtain a result in natural order (Assuming input is also in natural order). I also understood that this bit-reversal can be performed before or after the FFT or IFFT processing.

However, when interverting the order of operations, from :

bit_reversal();
IFFT();

to :

IFFT();
bit_reversal();

I don't get the same results. I would expect some difference due to rounding (I am working in fixed-precision here). But I get really completely different results.

Unfortunately, I am not allowed to share the code here, so my question is not about the implementation of the bit-reversal or IFFT, but more about fundamentals that I may have missed : can we operate bit-reversal equivalently before or after the radix-2 IFFT, without changing the radix-2 algorithm ?

Thanks in advance.

To answer my own question : no you can't use the exact same algorithm.

In the case here of a decimation-in-time, radix-2 Cooley Tuckey, in-place FFT/IFFT algorithm, we have the choice between two implementations which are exactly equivalent :

The above uses bit-reversed inputs and produce natural order results.

The above uses natural order input and produces bit-reversed outputs.

You can see that the operations are the same : for each butterfly operators, same input data, same twiddle factors. It's just the address of that data (here, the number of lines from the top) that changes.

The same applies for IFFT. When applying the bit-reversal after the IFFT, I would have needed to change the IFFT algo so that it changes the addressing scheme to match the second algo.

To go a bit further on the topic : "ALGEBRAIC FORMULATION OF THE
FAST FOURIER TRANSFORM" from Tran-Thong. (http://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=6323745)


Mismatched offset when computing cross-correlation using FFT In MATLAB if I compute the cross-correlation of 2 arbitrary signals using the xcorr function I get exactly the answer I expect. If I then calculate the cross-correlation using element-wise fft multiplication, the resulting signal is incorrect (offset) by 1 element. Can anyone tell me what I'm missing? Here is my MATLAB script:

a = [0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0];
b = [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0];
elems = length(a);

% Calculate the cross-correlation coefficients using built-in function
[ac_coeffs, lags] = xcorr(a, b);
subplot(2, 1, 1);
plot(lags, ac_coeffs);

% Now do it manually by fft multiplication
a = [a, zeros(1, elems - 1)];
b = [b, zeros(1, elems - 1)];
af = fft(a);
bf = fft(b);
mult = af.* conj(bf);
result = ifft(mult);
subplot(2,1,2);
plot(result);

You're missing an fftshift to get the correct alignment:

plot(lags,fftshift(result));


Improving the separation of a signal into source and filter I am implementing a voice morphing system, but I am having difficulties getting a good quality output signal. The method I am using is LPC analysis and synthesis. In theory, the LPC process should split the signal into a source/filter representation by calculating a residual and a set of filter coefficients.

However, in practice, I am finding that the residual contains data which I believe should be contained in the filter.

What techniques are available for better separating an input signal into source and filter?

Cheers,
Ben

Have you considered using a higher order filter for your LPC analysis?  LPC analysis is known to be sensitive to the presence of noise, say additive, so it may be useful to preprocess/denoise the signal to have the best/cleanest signal going in. Are you using a microphone on a mobile device?

If those don't work, you should look into research on spectral envelope extraction, or alternative algorithms for voice-morphing.


relationship of frequency with range in Underwater acoustic channel I am interested to know what range of frequencies can be used/ or are effective for wireless communication, within the underwater acoustic channel. I do know a table exists from \cite{akyildiz2004challenges, akyildiz2005underwater, kheirabadi2013greedy}, which is reproduced below. I want to add another column in this Table to add the range of frequencies that would be effective for the given range of propagation distances. But cannot get an idea of how to get it done.

What determines the effective range of frequencies for a given frequency and distance in underwater acoustic channels are the frequency and distance dependent path loss, and the power spectral density (PSD) of the noise.

The path loss is given by

$$A(f,d)=d^{\kappa}\left[\alpha(f)\right]^d$$

where $d$ is the distance between the transmitter and receiver in Km, $f$ is the frequency in KHz, $\kappa$ is the spreading factor which is usually assumed 1.5 for practical consideration, and $\alpha(f)$ is the absorption coefficient which usually calculated using the Thorp's formula as (measured in dB/Km)

$$10\log_{10}\alpha(f) = 0.11\frac{f^2}{1+f^2}+44\frac{f^2}{4400+f^2}+2.75\times 10^{-4}f^2+0.003$$.

The PSD of the noise is the sum of the PSD of the noise sources in UWA channels, which are: turbulence, thermal, shipping, and wave noise. You can find the equation of each noise in the literature.

Having the path loss and the PSD of the total noise, you can find the 3-dB bandwidth as following: the narrow-band signal-to-noise ratio (SNR) is given by

$$\gamma(f,d)=\frac{P/A(f,d)}{N(f)\Delta f}$$

where $N(f)$ is the PSD f the total noise, $P$ is the transmit power, and $\Delta f$ is the frequency over which the channel is frequency-flat. For a given distance $d$, there is an optimum $f$ that maximizes the SNR (which basically minimizes the term $\frac{1}{A(f,d)N(f)}$). Let that frequency denoted by $f_d$. Then you find the range of frequencies for a given distance such that

$$\gamma(f,d) \geq \frac{1}{2}\gamma(f_d,d)$$.

This way you can find the range of frequencies that can can be used for UWA communication for a given distance.

See this reference for more information.


Chebyshev pole calculation I am designing a Chebyshev filter and I am required to calculate the poles of it. It requires $1\textrm{ dB}$ pass-band ripple, a ripple factor of $0.5088$ and an order of $3$.

The equations for calculating the real and imaginary poles are shown below:

\begin{align}
s_k &= \sigma_k + j\omega_k\\
\sigma_k &= \sin\left(\frac{(2k-1)\pi}{2n}\right)\sinh(v)\\
\omega_k &= \cos\left(\frac{(2k-1)\pi}{2n}\right)\cosh(v)
\end{align}

However, when I calculate the imaginary part I get the following values:

\begin{align}
\sigma_1 &= -0.96\\
\sigma_2 &= 0\\
\sigma_3 &= +0.96
\end{align}

However, according to my solutions I should get $+0.96$ and $-0.96$. This checks out with the values from this link.

Is there some sort of rule to invert the values?

As I've mentioned in a comment, the poles are not "ordered" in any way. You can come up with different versions of the formula for the pole locations, which will give you the same poles in different orders.

But anyway, looking at the formulae in your question, the first thing to notice is that there should be a negative sign in the formula for the real parts $\sigma_k$. Second, the values that you got must be the imaginary parts $\omega_k$ (and not $\sigma_k$ as written in your question). If I evaluate the formula for $\omega_k$ given in your question for $k=1,2,3$ I get the following values (for $n=3$ and $1$ dB pass band ripple):
$$\begin{align}\omega_1&=0.966i\\
\omega_2&=0\\
\omega_3&=-0.966i\end{align}$$

So I get a positive sign for $k=1$ and a negative sign for $k=3$.


Raised-Cosine Filter - calculating phase response I want to compute phase response of an raised-cosine FIR filter. I have trouble finding the complex transfer function of the filter in literature (there is a real valued transfer function of frequency).

I have found formula for magnitude for raised-cosine filter (paper “Raised Cosine Equalization Utilizing Log Scale Filter Synthesis”, D. McGrath, J. Baird, B. Jackson) but there is no mention of phase response. I have read the IEEE paper “Generalized Raised-Cosine Filters” by N. S. Alagha, P. Kabal, but did not find the clear way how to calculate the phase response. In the first mentioned paper the phase is not linear and it's said that it is minimum phase filter. Here is the plot from McGrath at al. paper, page 8:

What is the way to calculate phase response for an raised-cosine FIR filter? I don't want to use MATLAB, I need to implement the solution in C++.

Just like the sinc pulse, a raised cosine pulse with zero delay (centered around $t=0$) has a real Fourier transform. You can find its expression in section 11.4 of "Telecommunications Breakdown" by Johnson and Sethares (a free draft version is available here).

For a time-shifted pulse, you can apply the time-shift property of the Fourier Transform.


Calculation frequency response of digital filter with known structure Short question
 What are main stages (steps) of calculation frequency response of digital filter by their structure?

Detailed question
Let suppose that there is discrete FIR filter with known structure which is implemented in programming language  (for instance with structure which shown in picture):

              +-----+   +-----+
              |     |   |     |
x (k) >---+---|  T  |---|  T  |
          |   |     |   |     |
          |   +-----+   +-----+
          |      |         |
          |   +-----+      |
          |   |     |      |
          |   | x 2 |      |
          |   |     |      |
          |   +-----+      |
          |      |         |
          |   +----------------+   +-----+
          |   |       ___      |   |   1 |
          \---|       \        |---| x - |---> y (k)
              |       /__      |   |   4 |
              +----------------+   +-----+

It is possible to pass into the filter input signal (as vector of integers which describe magnitude), like this:

in [0, 0, 0, 0, 4, 0, 0, 0, 0, 0]

Which describe signal with following magnitude:

             |
x (k)        |
             |
   . . . . . | . . . .

k  9 8 7 6 5 4 3 2 1 0

And get corresponding output signal, like this:

out [0, 0, 0, 0, 1, 2, 1, 0, 0, 0]

Which describe signal with following magnitude:

y (k)
           |
   . . . | | | . . . .
k  9 8 7 6 5 4 3 2 1 0

Question is how to calculate frequency response of the filter (by stages)?

Notes
With known input and output signals we can calculate transfer function of the filter by their images of Laplace transform (z-transform for discrete signals):
1. $$X(z) = \mathcal{L} \{x(k)\}$$
2. $$Y(z) = \mathcal{L} \{y(k)\}$$
3. $$H(z) = Y(z) / X(z)$$
And after if we get Fourier transform from transfer function we will get AFC:
4. $$FR = \mathcal{F} \{H(z)\}$$
Is it correct?

For the given system you can write down the input-output relation as

$$y[k]=\frac14\left(x[k]+2x[k-1]+x[k-2]\right)\tag{1}$$

because $T$ (or $z^{-1}$) denotes a delay element, which delays its input by one sample interval.
The $\mathcal{Z}$-transform of (1) is (assuming zero initial conditions)

$$Y(z)=\frac14\left(X(z)+2X(z)z^{-1}+X(z)z^{-2}\right)=\frac{X(z)}{4}\left(1+2z^{-1}+z^{-2}\right)\tag{2}$$

From (2) you get the system's transfer function

$$H(z)=\frac{Y(z)}{X(z)}=\frac14\left(1+2z^{-1}+z^{-2}\right)\tag{3}$$

Since the system is stable (any FIR filter is), the frequency response can be obtained by evaluating the transfer function on the unit circle $z=e^{j\omega}$:

$$H(e^{j\omega})=\frac14\left(1+2e^{-j\omega}+e^{-2j\omega}\right)\tag{4}$$

You also could have arrived at (4) by writing down the impulse response from (1)

$$h[k]=\frac14\left(\delta[k]+2\delta[k-1]+\delta[k-2]\right)\tag{5}$$

and taking the (discrete-time) Fourier transform, which also results in (4).

Finally, a few notes concerning misconceptions in your question:

$\mathcal{L}$ usually denotes the Laplace transform, which is defined for continuous functions. $X(z)$ is the $\mathcal{Z}$-transform of the sequence $x[k]$, which can be written as $X(z)=\mathcal{Z}\{x[k]\}$.
The frequency response of a linear time-invariant (LTI) system is the Fourier transform of the system's impulse response: $H(e^{j\omega})=\mathcal{F}\{h[k]\}$, not the Fourier transform of $H(z)$ (which you seem to think, judging from the last equation in your question).


Band Limit in Gaussian Pyramid In the paper The Laplacian Pyramid as a Compact Image Code,

it says The pyramid algorithm reduces the filter band limit by an octave from level to level.
How is the conclusion drawn? Is this related to weighting function and downsampling width and height by 2? How does weighting function and downsampling operation affect the filter band limit?

Edit:
From one document, it shows the frequency response of binomial filters

It's due to the weighting function used in the Gaussian pyramid downsampling.

Lets take a 1D patch(1D patch because it makes our calculation easy and also the kernel is separable) with p10,p11,p12,p13,p14,p15 in the Level 1.

Assume maximum gradient possible is Max - Min between any two pixels.

The pixel p20 in level 2 is formed by influence of pixels p10,p11,p12,p13,p14 in level1 and pixel p21 if formed by influence of pixels p12,p13,p14,p15,p16. The corresponding image(I'm not not able to upload image).

The weights of the kernel  are

w(0) = a,

w(-1) = w(1) = 1/4

w(-2) = w(2) = 1/4 -a/2.

  For Gaussian pyramid a = 0.4(approximately)

  w(0) = 0.4

  w(-1) = w(1) = 0.25

  w(-2) = w(2) = 0.05

The equation for new pixel p20

  p21 = w(-2)*p10 +w(-1)*p11+w(0)*p12 +w(1)*p13+w(2)*p14

Substituting the values for the weights we get,

  p20 = 0.05*p10 + 0.25*p11 + 0.4*p12 + 0.25*p13 + 0.05*p14

Similarly value of pixel p21

  p21 = 0.05*p12 + 0.25*p13 +0.4*p14 +0.25*p15 + 0.05*p16

Now the difference between pixels p20 and p21 is,

  p20 - p21 = 0.05*p10 + 0.25*p11 +0.35*p12 -0.35*p14 -0.25*p15+0.05*p16

So for the difference to be maximum p10 = p11 =p12 = Max and p14=p15=p16 = Min

  p20 -p21 = 0.7(Max-Min)

  p20- p21 = (Max -Min)/√2

Since the filters are separable, applying the kernel in another direction over the image again will reduce the pixel difference by a ratio of √2.
The maximum gradient in x and y direction in next level of Gaussian pyramid is (Max - Min)/2

The maximum gradient possible is reduced by half or in frequency domain reduces the frequency by half(an octave).


Signal Power estimation I want to add a White Gaussian noise to noise-free signal given desired SNR value. To do that, I am looking for a way to estimate the power of the signal. The noise-free signal looks like the following:

Do you have any clue about how to estimate the power of the signal ( mean of the signal, the integral of the difference between the decreasing linear line  and the concave down decreasing line, ....)?
Many thanks in advance!!!

The plot has the power already measured for each sample. If the raw data samples (and not the power results shown here, then the average power can be direclty measured (given by the variance of the waveform). We may think that simply averaging the power measurements in the plot would give us the average signal power--- This would be the average of the squares which is not the same as a true power measurement given by the square of the average. It will be close but such an approach is highly dependent on the statistics of the waveform. For example, if the waveform was a complex Gaussian signal, the average of the squares would under-estimate the true power by -1.05 dB. This is detailed nicely in the app note by Agilent (formerly HP) # 1303 available at this link http://hpmemoryproject.org/an/pdf/an_1303.pdf .
Typically SNR would be determined from the noise power within the occupied bandwidth of the signal while a white noise contribution is spread evenly over the sampling bandwidth. If the signal is oversampled, then you would want to increase the white noise source accordingly to maintain a target spectral density within the bandwidth of the signal.


This is an expression for the computation of kurtosis. However, I don't understand what the subscript '4x' or the parameter (0,0) stand for. Could anyone explain ?

For a stationary process $x$, moments and cumulants, when they exist, depend on relative time lags, hence $n-1$ variables when the moment or cumulant is of order $n$ (see for instance Basic Signal Processing and Statistics, p. 446 sq.):

$ \mathcal{M} \{ x(k),x(k+\tau_1),x(k+\tau_2),x(k+\tau_3) \} = \mathrm{E} \{ x(k)x(k+\tau_1)x(k+\tau_2)x(k+\tau_3) \}\,. $

A standard notation is $m_{4x}(\tau_1,\tau_2,\tau_3)$, where $x$ denotes the  process, $4$ the $n$-th order, $\tau_1,\tau_2,\tau_3$ the lags, with a triple $0$ for the $0$-lag cumulants. The first $0$-lag cumulants are classicaly named variance, skewness and kurtosis. Your version seems a normalized kurtosis, which could be written, with the above notations:
$$ \kappa_{4x}(0,0,0) = \frac{m_{4x}(0,0,0)}{m_{2x}^2(0)}\,. $$
I am quite surprised thant you only have two zeroes in your parenthesis.

On the  origin side, kurtosis comes for a greek word meaning arched, curved, measuring the tailedness of a distribution.


I want to perform aec on a wav file and calculate erle I am trying to evaluate the performance of an AEC by comparing the i/p and output wav file. I want to calculate echo return loss enhancement to predict how much(db) echo has been attenuated but i am not sure what would be right approach for that or is it possible or not, can we calculate in python?

You could use cross correlation for this purpose. For an uncorrelated signal $x(t)$, $\textrm{corr}(x(t), x(t + \Delta t)) = c \delta(\Delta t)$, i.e. an delayed impulse. This can be expanded to comprise multiple echoed copies
$$ \textrm{corr}\left( \sum_{i=0}^{N} a_i x(t + \Delta t_i), \sum_{j=0}^{N} b_j x(t + \Delta t_j) \right)  =
c \sum_{i=0}^{N} \sum_{j=0}^{N} a_i b_j \delta(\Delta t_i - \Delta t_j) $$
I remember from highschool physics books that if a delay is less than a certain amount, associated with the low frequency human hear response, (say 0.1 second), then it is not echo, it is reverberation. More objectively we could say that the audio signal will have a cross correlation that vanishes for lags larger than $L_x$.
You can have a measurement of the echo energy by the sum of the cross correlation with lag longer than $N_e$ from the maximum correlation. The maximum correlation will also provide a reference for synchronization of the two audio signals.
Putting this together you end up with a code like this
def measure_echo(a1, a2 = None):
    corr = signal.correlate(a1, a2 if a2 is not None else a1)
    lags = signal.correlation_lags(len(sig), len(sig_noise))
    # find the peak correlation (synchronize the signals)
    center = np.argmax(np.abs(corr))
    return np.sum(np.abs(corr**2))/np.sum(np.abs(corr[center-Ne+1:center+Ne]))


FFT of a pure sine wave is not a single peak For fun, I computed the FFT of a pure sine wave. I chose chose the sample length to be an even multiple of the signal period, so that I don't see windowing-effects. Here is the Matlab code that I used:

% Signal properties
f_signal = 1000;

% Sampling properties
f_sample = 4000;
T_total = 100;

% Calculate signal vector
t = 0:1/f_sample:T_total;
N = length(t);
V = sin(2*pi*f_signal*t);

% Calculate FFT
f_spectrum = (0:N-1)/N*f_sample;
V_spectrum = abs(fft(V))/(N/2);

% Plot
figure(2)
semilogy(f_spectrum, V_spectrum)
set(gca, 'XLim', [f_signal-5, f_signal+5])
xlabel('Frequency in Hz')
ylabel('Normalized amplitude')

But instead of a single peak, I get the following result:

This strikes me as odd for two reasons:

The effect is really huge. The center peak has an amplitude of 0.9 instead of 1, and the next bin to the right has an amplitude of 0.3 instead of 0.
While I would expect that limited numerical precision broadens my peak a bit, I am quite surprised that the effect is so silky smooth. I always thought of rounding errors as "noise" of the least significant digit. Noise usually leads to jarry spectra.

I observed that if I choose a slightly lower frequency (999.9975) then the effect is minimal. What is going on here?

Just as I wanted to click the "ask" button, I realized that I had run into a fencepost problem.

My time vector is too long by 1 point. This means that the first and last point of my voltage vector are the same, which in turn means that my periodic signal has a phase shift as it repeats. This also explains why the phenomenon went (almost) away as I decreased the frequency.

If I instead define my time vector properly as

t = 0:1/f_sample:T_total-1/f_sample;

then I get the result that I expect: A single peak and a tiny bit of noise around it.


Combining two wave generation formulas does not work as expected (MATLAB/OCTAVE) I'm trying to convert sawtooth wave calculation formula to Octave/Matlab language. Here are the original formulas. My Desmos take is here. I've got formulas for triangle wave and square wave working but, sawtooth won't work in either way like showed in Desmos sheet (two last formulas there).
Here's my Octave code (active saw calculation is simplified version of the above equation in comments):
function saw = saw_wave_trig(fs, f0, s, m, w)
% fs = sample rate
% f0 = fundamental frequency
% s  = 'sharpness'
% w  = 'width'
% m  = 'tooth direction'
% example: saw = saw_wave_trig(44100, 440, 0.03, -2, 1);

t = 0:1/fs:1;
##saw = m .* ((1.+ triangle_wave_trig(fs, (2*(f0/w)+1)/4, s, w) .* ...
##              square_wave_trig(fs, f0/2, s, w))./2) ...
##     .- sign(m);

##saw =  m .* ...
##            (1 .+ ...
##                  ((1.-(2.*acos((1-s) .* sin((2*pi*((2*(f0/w)+1)/4)).*t')))./pi) .* ...
##                       (2 .* atan(sin((2*pi*(f0/w/2)) .* t')./s))./pi)) ./ ...
##                         2 .- sign(m);

saw = (4.*asin((s-1).*sin((pi*t'.*(2*f0+w))./ ...
        (2*w))).*acot(s*csc((pi*f0.*t')./w)))./pi^2-sign(m)-1;

audiowrite('test_saw_trig_wave.wav',saw, fs);
plot(t, saw);

which (the commented calculations as well) produces kind of morph from negative side triangle wave to correct looking sawtooth wave, as seen in attached plot (Audacity used):

So, as my knowledge in math and Octave/Matlab is quite limited, there sure is some basics done wrong in my Octave/Matlab code. Maybe the multiplication of triangle and square needs to be done differently?
Any help in solving this is appreciated.
NOTE! This method is part of the experiment I'm doing with various methods out there to synthesize sawtooth (and other type) waves for audio use.
EDIT:
Here are those triangle_wave_trig() and square_wave_trig() functions (f1() and f2() in Desmos sheet) which of both works correctly:
function tri = triangle_wave_trig(fs, f0, s, w)
% fs = sample rate
% f0 = fundamental frequency
% s  = 'sharpness'
% w  = 'width'
% example: tri = triangle_wave_trig(44100, 440, 0.03, 1);

  t = 0:1/fs:1;
  tri = 1 .- (2 .* acos((1 - s) .* sin((2 * pi * (f0 / w)) .* t')) ./ pi);

  audiowrite('test_triangle_wave.wav',tri, fs);
  plot(t, tri);

function sqw = square_wave_trig(fs, f0, s, w)
% fs = sample rate
% f0 = fundamental frequency
% s  = 'sharpness'
% w  = 'width'
% example: sqw = square_wave_trig(44100, 440, 0.03, 1);

  t = 0:1/fs:1;
  sqw = 2.0 .* (atan(sin((2 * pi * (f0 / w)) .* t') ./ s) ./ pi);

  audiowrite('test_square_wave_trig.wav',sqw, fs);
  plot(t, sqw);
  hold on;
  plot(t, square(2*pi*f0*t'));

sqw = square_wave_trig(44100, 440, 0.03, 1); looks like this:

and generated wave file looks OK (same with the triangle wave function).

OK, got this solved by @omersayli (dsprelated).
I had put the variable t to a wrong location in both (triangle and square wave) formulas. Here's the correct formula for calculation of saw:
saw = -(2.*m.*asin((s-1)*cos((pi*f0*t)/w)).*acot(s.*csc((pi*f0*t)/w)))./pi^2.-sign(m).+m./2;


Techniques to reduce delay in detecting face from video using vision.CascadeObjectDetector in matlab I wrote matlab code for face detection.In my code it is detecting face for first 100 frames and it crop faces from each frame and saves it in database folder.Problems iam facing

1.Detecting face frame by frame is very slow.Is there any idea to run faster since i have to work on 4000 frames.

2.In my database folder it has to show 1 to 100 face images but it is not showing 11th and 12th face images directly it showing 13th face image after 10th image.23rd face image is blurr.Likewise so many images are missing and some are blurr.Last image number it is showing as 216.But total 106 face images are there in database folder.In that 12 images are blurr.Remaining are correct images.

clc;
  clear all;

  obj=vision.VideoFileReader('basu.avi');

  for k=0:99

  videoFrame      = step(obj);
  %using viola-jones algorithm
  FaceDetect = vision.CascadeObjectDetector;

 %FaceDetect
  BB = step(FaceDetect,videoFrame);
  %BB
   figure(2),imshow(videoFrame);

  for i = 1:size(BB,1)
  rectangle('Position',BB(i,:),'LineWidth',3,'LineStyle','-','EdgeColor','r');
  end

  %crop faces and convert it to gray
  for i = 1:size(BB,1)
  J= imcrop(videoFrame,BB(i,:));
  I=rgb2gray(imresize(J,[292,376]));

  %save cropped faces in database folder
  filename = ['G:\matlab_installed\bin\database\' num2str(i+k*(size(BB,1))) '.jpg'];
  imwrite(I,filename);
  end
  end

There are a few of things you can try:

Definitely move FaceDetect = vision.CascadeObjectDetector; outside of
the loop. You only need to create the face detector object once.
Re-creating it for every frame is definitely your performance
bottleneck.
vision.VideoFileReader returns a frame of class 'single' by default.
If you change the output data type to 'uint8', that    should speed
up    the face detector.Use obj=vision.VideoFileReader('basu.avi',    'VideoOutputDataType',
'uint8');
vision.VideoFileReader can also do the conversion to grayscale for
you. Use obj=vision.VideoFileReader('basu.avi',
'VideoOutputDataType', 'uint8', 'ImageColorSpace', 'Intensity'); This
may be faster than calling rgb2gray.
Try limiting the size of the faces being detected using 'MinSize' and
'MaxSize' options of vision.CascadeObjectDetector    and/or try
downsampling the frame before detecting faces.


Using Pilot tones to estimate carrier frequency offset in ofdm I have been studying carrier frequency offset estimation techniques in OFDM systems. I have understood several of them. when using Pilot tones, as stated in paper (1), the process involves two methods namely: acquisition and tracking modes. a large range CFO is estimated in the acquisition mode while in  the tracking mode, only fine CFO is estimated.
I have confusion, What are the nature of those Pilot tones? how are they being generated?
the authors of the paper wrote equation (18) without any derivation? can I get the system model for this technique?
the questions may sound easy, but I really need it. Thanks in advance.

Reference:
(1) Praween Kumar Nishad,  P. Singh, "Carrier Frequency Offset Estimation in OFDM Systems", National Institute of technology
http://dspace.nitrkl.ac.in/dspace/bitstream/2080/1928/1/ID-451_praween.pdf

I am familiar with 802.11a, so I'll answer your questions with regard to how 802.11a does it.

What are the nature of those Pilot tones?

There are 64 sub-carriers in 802.11a.  Four of them are used as the pilot channels.  The pilots themselves are simply known data- if I remember correctly they are either (1 + i) or -(1 + i).

how are they being generated?

The transmitter just puts the appropriate data into the pilot sub-channels and then ifft's the entire symbol (i.e. with all of the other sub-channels).

the authors of the paper wrote equation (18) without any derivation?

This isn't really a question.

can I get the system model for this technique?

What do you mean by a system model?  We don't normally provide code.


Laplacian Filter Implementation in MATLAB I implemented a Laplacian filter for the Lena image, but I get an unexpected output.

For input:

I get output:

I used the mask

M =

     0    -1     0
    -1     5    -1
     0    -1     0

I is my Lenna Image rows and columns are I.rows and I.columns.

And this is my MATLAB function:

function IM = customMaskImage(I, rows, columns, M)
IM = zeros(rows,columns);
[mask_rows, mask_columns] = size(M);

for i = 1+(mask_rows-1)/2:(rows - (mask_rows-1)/2)
    for j = 1+(mask_columns-1)/2:(columns - (mask_columns-1)/2)
        tmpMat = I(i-(mask_rows-1)/2:i+(mask_rows-1)/2,j-(mask_columns-1)/2:j+(mask_columns-1)/2);
        IM(i,j) = sum(sum(customMatConv(tmpMat, 3, 3, M)));
    end
end

end

customMatConv.m is:

function IM = customMatConv(I, rows, columns, M)
IM = zeros(rows,columns);

for i = 1:rows
    for j = 1:columns
        IM(i,j) = I(i,j) * M(i,j);
    end
end

end

What am I doing wrong?

You should use conv2() or imfilter() with your filter.
Unless you want to implement it by hand,


Why does AFSK only trasmit in binary? According to Wikipedia, audio frequency-shift keying (AFSK) is done by using binary FSK, i.e. 2 tones of different frequencies, a tone for the 1 bit and a tone for the 0 bit.

But isn't it possible to choose more than 2 frequencies (any frequency between 0 and half the sampling rate) to transmit more than binary data?
Why don't we use say 10 tones of different frequencies, to transmit in decimal, or say 256 tones to transmit in chars?  Wouldn't that make transmission much more efficient?

The three main reasons appear to be these:

1) In AFSK, it's not just jumping back and forth between frequencies. The tones must also be continuous phase. In other words, when the tone changes, there can't be any jump in phase. For example, if you're sending a 1200 hz tone and the waveform is at its peak when you switch to 2200 hz, the waveform is still at its peak, so it can't be just any two frequencies. It has to be at the peak. It can't start back at zero, or halfway through; it can't start anywhere but at the peak. Additionally, the length of each tone must be the same. Because of these two points, the amount of wavelengths that share compatible periods is finite. For instance, if you use 1000 hz, your next compatible is something like 1500 hz (my apologies if my numbers are off; I'm worn out).

2) As a follow-on from 1, since AFSK is audio, the range of compatible frequencies is not that high.

3) The amount of bandwidth required for the transmission increases for each additional tone, making data efficiency greater, but power and bandwidth efficiency less.

That doesn't meant that multiple frequencies can't be used; in fact MFSK has been used since the mid-60s. However, MFSK has several drawbacks that make it untenable in some situations. Delay spread, when broadcasting to multiple locations, and fading make it less desirable for some applications.


Interpreting this plot cross-correlation I have an input signal $x(t)$ that is a white Gaussian random signal with mean 0 and variance 1. The signal $y(t)$ is the output of a linear filter with impulse response $$h(t) = |\operatorname{sinc}(t − 2)|.$$ I have plotted the cross-correlation and auto-correlation of the two signals, but have a hard time interpreting what the plots mean. For the cross-correlation, what do all the different spikes mean? What does the peak near 0 mean?
For the auto-correlation of $x$ as well, why is there noise but a peak at 0?
What about the auto-correlation of $y$?

Input and output signals:

Cross-correlation:

Auto-correlation:

Autocorrelation of $x$: at zero shift, the two signals being correlated ($x(t)$ and $x(t)$) are identical, hence a large correlation. At any other shift ($x(t)$ and $x(t+\tau)$), the correlation is near zero (no correlation). Some noise there is expected.

Autocorrelation of $y$: at zero shift, as always, there is maximum correlation of the two signals. Because the signal is smooth, correlation changes slowly with increasing shift. At various shifts peaks line up an the correlation has a local maximum.

Crosscorrelation of $x$ and $y$: you’re not showing the magnitude of the y-axis, I presume it’s small because the signals are not correlated. At some shift (it just happens to be close to zero, this is coincidence) there is the absolute maximum because it has to exist somewhere. Again, because $y$ is smooth, the correlation changes smoothly with the shift. But in general you’re looking at smoothed noise here.

If you computed the crosscorrelation by padding the input with zeros (rather than using the circular periodicity imposed by the Fourier transform) then you would expect the peak to be close to 0, as the correlation for larger shifts will include more of the zero values in the padding (the signals overlap less, hence the sum of their point-wise multiplication is closer to 0).


Entropy of a dynamical system and source code length The entropy of a source $H(S)$, gives you the average codeword length to encode a given source alphabet. i.e. it is the average number of bits per symbol required to encode the information in the source. If: $p_0=3/5$ and $p_1=2/5$ (close probabilities) with $H=0.971 \textrm{ bits/symbol}$, then there's more information than compared to the 2nd case where $p_0=3/4$ and $p_1=1/4$ yielding $H = 0.8113 \textrm{ bits/symbol}$ to encode. Following are the problems that I am facing and I will be very grateful for a detailed answer.

In general, if the source entropy is $H$ then what should be the minimum number of bits if the alphabet is $0,1$ ?
In case of discrete chaotic map say Logistic map, the entropy is $\log \lambda$ where $\lambda$ is the largest Lyapunov exponent. What should be the length of the time series in binary or what should be the length of the symbolic representation in binary of the initial condition?

I can try to answer your first question (but I'm not sure I understand exactly what you mean).

A source produces messages with average information $H$. You seem to have a binary source, with two messages $m_1$ and $m_2$.

Now you want to encode the messages to a sequence of 1s and 0s (subsequently called "bits"). This is called "source coding". Even though we know that, in principle, we should be able to encode $N$ messages using $NH$ bits (for large $N$), it is not always obvious how to do so. There are algorithms such as Huffman and Lempel-Ziv that can find good (even optimal) encodings. Note that all good coding algorithms encode messages in groups, not individually, and they approach $NH$ only when $N$ is large.

If your question is "what's the smallest number of bits that I need to encode $N$ messages from a source with entropy $H$", the answer is $NH$ (in theory), and somewhat more in practice, depending on the optimality of your encoding algorithm and the value of $N$.

Example 1: A source produces 8-bit messages $b_7b_6b_5b_4b_3b_2b_2b_1b_0$ uniformly at random, but $b_7$ is always 0. The entropy is 7 bits/message. Let's say we want to transmit 100 messages. Even though the source produced 800 bits, we only need 700 to transmit the messages.

Example 2: You flip a biased coin with $p(\text{heads})=0.3$ 1000 times. The entropy of the experiment is 0.88 bits/flip. To convey the result of the 1000 flips, you'll need on average 880 bits. Some particular results will need more bits, some will need less. To figure out how many bits you need for a particular sequence of flips you need to specify your encoder algorithm.

Example 3: You have a biased 3-sided die with probability distribution $p(1)=0.5$, $p(2)=0.25$, $p(3)=0.25$. The entropy is 1.5 bits/throw. It's pretty easy to come up with an encoder:

result  encoding
1    -> 0
2    -> 10
3    -> 11

You should be able to convince yourself that, on average, this encoding results in 1.5 bits per throw. However, let's say that you're unlucky and you get 1000 '3' in a row. That particular sequence will need 2,000 bits to transmit. That's why the entropy can only predict the average sequence length, not the length you'll get for a particular set of messages.


Conversion of transmitter power from Watts into dB I've been practicing problems relating to the dB conversion and came across this interesting example problem.

In the above example Tx power is defined in Watts and the conversion from the Watts to dB is done. However, so far according to my understanding of the definition of dB, is defined in relative to receiver power too. But in the above example the calculation was done just by considering the Tx power.

Can some one enlighten me? If I'm understanding it wrong or the example is wrong?

Thanks in advance.

Thanks to @MBaz who pointed out my mistake: dBm and dBW are unitless, but they relate to  a reference of 1 mwatt, and 1 Watt.

The question has an error, as well as the answer. Watts does not convert to dB, it converts to dBW. (i) should be as follows:

find (i) the transmit power in dBm and dBW.

The reason behind that is that dB is a ratio, while dBm and dBW are ratios with a reference to 1 mWatt and 1 Watt, respectively, and you can use them to recover an absolute value of the power.

For example, the signal to noise ratio is a measure of how strong a signal is compared to the noise. In the linear scale it is written as:

$\frac{\text{Signal Power (in watts or miliwatts)}}{\text{Noise Power (in watts or miliwatts)}} = SNR (\text{unitless})$

or in logarithmic domain:

$\text{Signal Power (in dBW or dBm)}-\text{Noise Power (in dBW or dBm)} = SNR (\text{in dB})$

So the correct answer to (i) is

TX power = $10\log(50)=17\text{ dBW}=17+30=47\text{ dBm}$


Power of a Discrete time signal I was learning on power of signal and i went through different formula in internet .

I get different results when i go through internet .

But in other site average power is  written as

I dont know if both are different formula for different condition .

Thanks in advance

The power of a discrete-time signal $x[n]$ is given by

$$P_x=\lim_{N\rightarrow\infty}\frac{1}{2N+1}\sum_{n=-N}^{N}|x[n]|^2$$

which is identical to the first formula in your question. The second formula is wrong and does not even make any sense due to the infinite limits of the sum.

I'd advise you to use text books or high-quality internet resources such as MIT OpenCourseWare.


Matlab: Correct method of generating Rayleigh fading channel A Rayleigh flat fading channel has 1 tap and so there will be no intersymbol interference.
if

N= 1000;
 Then,
h=(1/sqrt(2))*(randn(N,1)+1i*randn(N,1));

will generate 1000 coefficients which are complex Gaussian of zero mean.

Question 1: Is this a multi-path Rayleigh channel of number of taps = 1000?(Example taken from : http://www.mathworks.com/matlabcentral/fileexchange/30477-ber-curve-for-bpsk-in-rayleigh-channel/content/BER_BPSK_Rayleigh.m )

But, there are models such as Clarke's method of sinusoid and Jake's algorithm for modeling the Rayleigh channel as FIR. If the method in Question 1 and these algorithm create multi-path Rayleigh, then when do we use the Algorithm?

Question 2: In the code and explanation therein for Rayleigh flat fading code , a channel of length N is created but it is called flat. I am really confused as to what is the proper way to create multi-path Rayleigh fading and only flat Rayleigh fading (1 tap).

Please help.

Question 1:

Yes, it is a multipath Rayleigh Channel with $ 1000 $ taps. You need to know that when you generate complex coefficients, and the real part and the imaginary part of this coefficients are randomly chosen according to a Gaussian distribution, the magnitude will have a Rayleigh distribution and the phase will have a uniform distribution between $ 0 $ and $ 2 \pi $. The problem with this method, is that you don't have control on how the channel varies: the channel in a given instant of time will have no correlation with the channel in the previous instant. To add this correlation, we use the others methods as Clark's and Jake's models.

Question 2:

Look what he wrote in his website, in the assumptions section:

  The channel is flat fading – In simple terms, it means that the multipath channel has only one tap. So, the convolution operation reduces to a simple multiplication.

In his code, he is performing a multiplication, not a convolution, when applying the Rayleigh Channel.

He generates $ N $ coefficients, but what he is trying to do is simulate the channel variation in time, so each symbols is multiplied by an attenuation (inside the vector h). The vector h is not a vector of a filter coefficients, but a vector that stores the channel attenuation for each symbol.


Is it a correct way to calculate the Mean Square Displacement as function of time? I would like to know if I calculate the MSD correctly, here is a Jupyter snippet I am using:

import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
%matplotlib inline

x = np.random.randn(100)

df = pd.DataFrame(x, columns=['x'])
df['time'] = np.arange(0,len(df))
df['MSD'] = pd.Series((x - x[0])**2).rolling(window=10, min_periods=1).mean()

x_plt, = plt.plot(x, label='x(t)')
msd_plt, = plt.plot(df['MSD'], label='MSD(t)')
plt.legend(handles=[line_up, line_down])

Is there a smart way to average over the variances? The original definition of MSD has average inside, but I don't relay averaging here because I have only one particle.

According to the article on which you base your MSD computation, it's specified as a single number (an average over N items) which is defined as the mean (average) of the sum of squared-distances between positions of $N$ particles at time $t$ and given reference positions $x_n(0)$ for each particle at a reference time.

Assuming that you want to compute MSD over a 1D particle displacement of $N$ items, you would use the following Matlab/Octave code to achieve it:

N = 1024;                 % Number of particles
std1 = 2.34;              % set some variance
xpos = std1 *sqrt(3)*2*( randn(1,N) - 0.5); % N particles with positions.
xref = 0;                 % using same reference position for each particle.
MSD = (1/N)* sum ( (xpos-xref).^2 );  % compute MSD of N particles at time t

figure,plot(xpos);        % plot positions of each particle at time t
hold on;                  % plot MSD of those N particles at time t.
plot(MSD*ones(1,N));      % MSD*ones(1,N) is used to plot N points for MSD
title('MSD over particle positions at time t');

Assuming that you also want to plot MSD as a function of sampled time, you could then then use the following code to achieve it:

N = 1024;                 % Number of particles
M = 128;                  % Number of MSD computations sample time points
std1 = 2.34;              % Set some variance, this can change over time.
xpos = std1 *sqrt(3)*2*( randn(N,M) - 0.5); % N particles for each M time samples.
xref = 0;                 % using same reference position for all time and each poarticle.
MSD = (1/N)* sum ( (xpos-xref).^2 );  % now MSD is an 1 x M vector

figure,plot(MSD);        % plot particle positions and MSD