Is this statement contradicting with the band-limiting theorem? The following sentences are from the footnote of page 240 in R.C. Gonzalez's Digital Image Processing (3rd):

  An important special case is when a function that extends from -Inf to
  +Inf is band-limited and periodic. In this case, the function can be truncated and still be band-limited, provided that the truncation
  encompass- es exactly an integral number of periods. A single
  truncated period (and thus the function) can be repre- sented by a set
  of discrete samples satisfying the sampling theorem, taken over the
  truncated interval.

Isn't these contradicting to the band-limiting theorem that a bandlimited signal cannot be also timelimited?

there is some choice of words i would do differently.

yes truncating the infinite bandlimited function makes it time-limited and then band-unlimited.  but truncating it to an integer number of periods still allows one to fully describe it with discrete samples.

a periodic function can be described with a Fourier series.  the Fourier coefficients will virtually entirely describe the periodic function.

a bandlimited and periodic function can be completely described with a finite set of Fourier coefficients.  also a bandlimited and periodic function can be completely described by sampling a single period with a sufficient number of samples that are spaced closer in time than $\frac{1}{2B}$ where $B$ is the bandlimit.


Finding remaining zeros of linear phase FIR filter Question given is: if $h(n)$ is a linear phase causal FIR of order $10$ with real coefficients, find the remaining zeros of this filter if the zeros given are

$$q_{1,2} = -2 \pm 2j$$

$$q_{3} = -\frac12 + \frac{j}{2} $$

$$q_{4} = -1$$

I know for sure that $-\frac {1} 2 - \frac j 2 $ has to be another zero, but there are 5 more zeros unaccounted for. How am I supposed to find them?

For linear phase FIR filters, each zero at z = z0 will have a matching reciprocal zero at z = 1/zo. And for real-valued coefficients each zero at z = zo will have a matching conjugate zero at z = *zo. Thus for linear phase real-valued coefficients, when you place one zero on the z-plane you determine the location of the other three zeros.


Measuring acoustic phase difference between two receivers I am trying to measure the acoustic phase difference between two sensors 30cm apart each other. The sampling frequency is 20kHz, and my interest frequency is about 1700 Hz.
The below is how I proceeded.

Take FFT of T = 1 sec time domain signal.
In spectrum, there are several peaks including the interest frequency.
Bandpass filter around the interest frequency (e.g. 1700+-100 Hz)
Reconstruct(or IFFT) the signal in time domain.
Take Hilbert transform and get phase in time domain
Do the same thing for the other receiver.
Compare Hilbert transform phase between two signals in time domain

When I did this, the phase difference between two signals is not constant, but changes with time, which is out of my expectation. Shouldn't the phase difference be constant if they have the same frequency peak? Is there any issue with the above procedure?

Keep in mind the following calculation:

$$ \frac{1130~\text{feet}/\text{s}}{1700~\text{cycle}/\text{s}} \approx 0.66~\text{feet/cycle} \approx 20.26~\text{cm}/\text{cycle} $$

Your delay distance is longer than your wavelength.  With an FFT approach, you will only get the delay as a fraction of a cycle, so you'll have to account for that.  If you can't frame your FFT on a whole number of cycles, you will need to estimate the frequency and calculate the phase shift from there.  Then, using your frequency you will need to calculate the actual delay distance.

Generally, if it is the same signal, you will get much better results using autocorrelation in the time domain.  Your next best bet is to make sure your FFT is aligned on a whole number of cycles, then you only need to calculate one bin.  If you can't do that, in my blog articles, I have frequency formulas from two or three bins and an method for determining the phase very accurately using the nearest three bins.  (See my profile page for a link to my blog).


Why are there vertical lines in the power spectrum? I took an image, did an fft2 to it in Matlab and then squared the abs of the fourier transform to get the power spectrum.

Something I notice about the power spectrum however is that there are vertical lines in it like this:

This is the original image which the power spectrum was calculated from, for reference:

What are these lines and why do I see such phenomenon from the transformed image? Do the lines imply anything about the original image and the fourier transform?

Let's interpret what we are seeing: the spectrum seems to be made of the pointwise multiplication of 2 known Fourier patterns:

a first one is the usual messy image Fourier transform
a second pattern that modulates the image spectrum and that has a $sinc$ shape.

Now, let's see what the theory provides us:

multiplication in the spectral domain is a convolution in the spatial domain => we are observing the convolution of an image by a blur kernel
the inverse Fourier transform of a $sinc$ is a box => the blur kernel has an horizontal box shape (we can infer that from the orientation of the zero bands).

Finally, let's summarize: the image was blurred by a linear (horizontal) kernel.
This is a typical example of motion blur, where the motion is horizontal.
This is confirmed by looking at the original image.

And as a bonus, here is a paper where special shapes of motion blurs are designed, where the blur shapes do not have zeros that will disturb the deconvolution process:
http://web.media.mit.edu/~raskar/deblur/CodedExpousreLowres.pdf


What is the Fourier transform of a multiplied cosine signal with rect? I tried to find the Fourier transform for this signal, but I have some doubt about the result:

$$g(t) = \cos\left(\pi\Delta{f}t\right) \text{rect}\left(\frac{t-\frac T2}{T}\right)$$

My calculation is:

$$ \bigg \lvert G(f) \bigg \rvert = \frac{1}{\pi}\bigg(\frac{2f\sin(\pi Tf)\cos(\pi T\Delta f/2)-\Delta f\cos(\pi Tf)\sin(\pi T\Delta f /2)}{4f^2-\Delta f^2}\bigg)$$

Could you please tell me is there any coeff. of $T$ at amplitude of the final result?

I think the only point in answering such overly specific questions is to give a more general view of how to approach such basic problems. So this is what I'll try to do here.

First of all, one should be able to recognize that the given signal can be viewed as an amplitude modulated rectangular signal. One should also know that a rectangular function in one domain of the Fourier transform is a sinc-function in the other domain. So from a first glance we should be able to tell that the resulting spectrum is composed of two sinc-functions, one shifted to the positive and the other to the negative frequency of the cosine. Finally, it should be observed that the frequency of the cosine is $\Delta f/2$ (not $\Delta f$). So we have two sinc-functions centered at $\pm\Delta f/2$. I strongly believe that such a qualitative analysis should always be done before trying to solve the problem in more detail.

Another global check reveals that the given result cannot be correct because it is real-valued. The true result cannot be real-valued since the given time-domain signal is not symmetric.

[EDIT: after I wrote this answer, the formula in the question was changed from $G(f)=\ldots$ to $|G(f)|=\ldots$. Now it's of course OK for the result to be real-valued, but note that it can still become negative, which it shouldn't if it were a valid expression for the magnitude of the Fourier transform.]

After all these sanity checks we might be interested in knowing the exact expression of the transform. This should be possible even without consulting a Fourier transform table. Since from the qualitative analysis we already know what the spectrum looks like, we only need to find the Fourier transform of the rectangular pulse, and then shift it appropriately. The Fourier relation $\text{rect}\Longleftrightarrow\text{sinc}$ is a very basic transform pair, and I think any DSP student or engineer should be able to get it right, including all necessary constants. One easy way to remember the correct relation is to memorize that a centered rectangular pulse of width $W$ and height $1$ in one domain corresponds to $W$ times a $\text{sinc}$-function with argument $W\cdot\text{variable}$ in the other domain, where "variable" is either $t$ or $f$:

$$\begin{align}\text{rect}(t/T)&\Longleftrightarrow T\,\text{sinc}(Tf)\\\text{rect}(f/B)&\Longleftrightarrow B\,\text{sinc}(Bt)\end{align}\tag{1}$$

where we need to use the definition $\text{sinc}(x)=\sin(\pi x)/(\pi x)$ (instead of the other common definition $\text{sinc}(x)=\sin(x)/ x$). Note that the relation $(1)$ is totally symmetric with respect to time and frequency, which is one advantage of using $f$ here instead of $\omega=2\pi f$.

Knowing $(1)$, we can immediately write down the transform of the given rectangular pulse:

$$\text{rect}\left(\frac{t-T/2}{T}\right)\Longleftrightarrow T\,\text{sinc}(Tf)e^{-j\pi Tf}\tag{2}$$

where the phase factor in $(2)$ comes from the shift by $T/2$.

From $(2)$ the final result is easily obtained by remembering the shifting property of the cosine (and the factor $\frac12$ that it introduces):

$$G(f)=\frac{Te^{-j\pi Tf}}{2}\left[\text{sinc}\left(T\left(f-\frac{\Delta f}{2}\right)\right)e^{j\pi T\Delta f/2}+\ldots\\\text{sinc}\left(T\left(f+\frac{\Delta f}{2}\right)\right)e^{-j\pi T\Delta f/2}\right]\tag{3}$$

I leave any final cosmetic operations on this result up to you.


Absolute and Square Summability & Integration of sinc function Is the sinc function both absolutely summable (L1 norm for Continuous time signals and l1 norm for Discrete time signals) and square summable (L2 norm for Continuous time signals and l2 norm for Discrete time signals) ?

Can anyone show the integration and summation calculations that calculates above norms ?

Neither the integral

$$I_1=\int_{0}^{\infty}\left|\frac{\sin(x)}{x}\right|dx\tag{1}$$

nor the sum

$$S_1=\sum_{n=1}^{\infty}\left|\frac{\sin(n)}{n}\right|\tag{2}$$

converge. So the sinc function (sequence) is not in $L_1$ ($l_1$). A proof of the fact that $I_1$ diverges can be found here. For the sum $S_1$, note that since $|\sin(n)|\ge \sin^2(n)$ (for $n\in\mathbb{R}$) we have

$$S_1\ge \sum_{n=1}^{\infty}\frac{\sin^2(n)}{n}\tag{3}$$

Here it is shown that the series on the right-hand side of $(3)$ diverges, so $S_1$ must also diverge.

The integral

$$I_2=\int_{0}^{\infty}\left(\frac{\sin(x)}{x}\right)^2dx=\frac{\pi}{2}\tag{4}$$

can be computed in many different ways. Two possibilities are shown in the answers to this question.

Also note that

$$\int_{0}^{\infty}\left(\frac{\sin(x)}{x}\right)^2dx=\int_{0}^{\infty}\frac{\sin(x)}{x}dx\tag{5}$$

The computation of the sum

$$S_2=\sum_{n=1}^{\infty}\left(\frac{\sin(n)}{n}\right)^2=\frac{\pi-1}{2}\tag{6}$$

is discussed in this question and its answers. Note that a relation similar to $(5)$ holds:

$$\sum_{n=1}^{\infty}\left(\frac{\sin(n)}{n}\right)^2=\sum_{n=1}^{\infty}\frac{\sin(n)}{n}\tag{7}$$

So clearly, the sinc function (sequence) is in $L_2$ ($l_2$).


Enhance image in matlab - exp operation - imageJ software-LPR I use imageJ software, imageJ has a EXP command (process >math >exp ).

original image:

output image processed by EXP command in imageJ software:

in this link explain how do that ImageJ: http://imagej.1557.x6.nabble.com/Process-gt-Math-gt-Exp-td5009733.html

How can I do that in matlab?

It is the matlab code which I wrote,Please correct my mistakes.

clear all;
close all;
clc;
img=imread('http://s6.uplod.ir/i/00594/fdg3933aznwz.jpg');
im=rgb2gray(img);
im=im2double(im);

figure,imshow(im)

    for i=1:size(im,1)
        for j=1:size(im,2)

          in = (im(i,j)+1);
          result(i,j) = 255*exp(in/46.01);

        end
    end
figure,imshow(mat2gray(result))

You have to learn to debug your programs. Try formula for some data:

>> 255*exp(255/46.01)
  ans =   6.5091e+004

So formula has bug (explanation in link to forum is valid. Read it carefully). Valid formula:

>> exp(255/46.01)
 ans =  255.26

But try calculation with real data from your image:

>> in = (im(10,10)+1)
in = 254
>> exp(in/46.01)
ans =  403.43    % very strange result!!!!
>> in/46.01
ans = 6          %what is it????

Reason is quite simple - type of in is uint8!!!

Valid formula:

result(i,j) = exp(double(in)/46.01);

result:

imshow(uint8(result));


About the order of a difference equation The question is related to a statement in a dsp book.  I am attaching the image of the statement.  In the attached image, is the equation 1.5 a difference equation of order 1 as stated in the book?  I think it is of order 0 since there is no feedback term i.e. no y(n-1).

To my understanding, the order of a difference equation is equal to the maximum delay, regardless if its delay on the input or output side. So, you have two first-order difference equations here.


explanation of Hybrid systems? I want to study Detailed explanation of hybrid systems?Which incorporate both continuous and discrete time signals & systems?
for example
In which a continuous-time input signal is transformed into a discrete-time output signal and vice versa. An analog to digital (A=D) converter is an example of a hybrid system.
Where can i find detailed explanation to quench my curiosity/thirst of this topic , i have tried searching google but no useful result

If you're not interested in specific practical aspects of A/D conversion, but if you want to learn basic theory concerning sampling and digital (discrete-time) processing of analog (continuous-time) signals, I'd recommend that you read and study the chapter Sampling of Continuous-Time Signals in Oppenheim en Schafer's book Discrete-Time Signal Processing.


Feedback systems & oscillations The transfer function of feedback system is:
$$
\frac{V_{\rm out}}{V_{\rm in}} = \frac{A}{1+Af}
$$
Where $A$ is the open loop gain, and $f$ is the feedback gain.

Now for oscillation to happen, $Af$ term (loop transmission) must be $-1$, which means the system is configured to be on a pole, because the denominator goes to zero. I have got couple of questions on this one.

For oscillation to happen (steady state), the pole must be on $y$-axis (frequency axis) in the $s$-plane. I believe this is an assumption in the oscillation theory not reflected by the above condition that loop transmission is $-1$. Because loop transmission is $-1$ in all poles in the system. So it is up to the designer to make sure that the circuit operates on a pole on $y$-axis. Is this correct?
If I ignore all the Laplace transform thing, and just look at the transfer function
$$\frac{V_{\rm out}}{V_{\rm in}} = \frac{A}{1+Af}$$
On the pole, the output is supposed to be very high (infinity). If I only consider this equation on standalone basis, then how come we have a steady state oscillation, the output voltage is supposed to be very high as denominator is zero. My explanation is that it does not happen, because the supply voltage to the Op-Amp/MOSFET (or whatever we use to implement the circuit) restricts that output voltage. Is this correct?

I think that the problem is that you are considering that $$ T = \frac{A}{1+Af} $$ is the actual transfer function of your system (an analog amplifier I guess) and this can be miss leading because that expression is just and approximation for the region of frequencys were the amplifier response is flat. If you want to make assumptions regarding the stability of the system you should work with the full Laplace transform, this is: $$ T(s) = \frac{A(s)}{1+A(s)f(s)} $$
In order to oscillate at a given frequency, the system needs to have conjugate poles on the imaginary-axis at that frequency, but in a practical case, if the system has poles on the right half-plane, as you said, the supply voltage of the amplifiers will move this poles into the j-axis

I don't have enough reputation to comment, but I think that this question that I asked a few weeks ago may be helpful.
https://electronics.stackexchange.com/questions/279365/positive-feedback-with-constant-transfer-functions


Windows Audio WAVEFORMAT buffer BYTE to 16-bit int I am trying to understand audio-buffer processing at the low-level for the first time.

I have the following C++ code copying WAVE-format byte data from capture buffer of Window Core Audio (*ppData in that documentation), into an "AudioSink."

Could someone explain what's happening here (the "WHAT'S GOING ON HERE" in the code)?
I'm confused as to why this works, how the number of channels, number of frames, and the conversion between BYTE and int16 factor into the calculation?  I need to eventually convert this into float data so I need to understand what is happening!

Thank you!!!

int AudioSink::CopyData(const BYTE* Data, const int NumFramesAvailable)
{
    // Data is the ppData returned pointer from IAudioCaptureClient::GetBuffer

    std::lock_guard<std::mutex> lock(m_mutex);

    AudioChunk chunk;
    if (Data == NULL)
    {
        chunk.size = 0;
        m_queue.push(chunk);
        return 0;
    }

    // *****  WHAT'S GOING ON HERE  *******
    chunk.size      = NumFramesAvailable * 2;
    chunk.chunk     = new int16[chunk.size];
    int multiplier = sizeof(int16) / sizeof(unsigned char);
    std::memcpy(chunk.chunk, Data, chunk.size * multiplier);
    bool nonZero    = false;

    for (int i = 0; i < chunk.size; i++)
    {
        if (chunk.chunk[i] == -1 || chunk.chunk[i] == 1)
        {
            chunk.chunk[i] = 0;
        }
        if (chunk.chunk[i] != 0)
        {
            nonZero = true;
            UE_LOG(LogTemp, Log, TEXT("NumFramesAvailable: %d, Sample number %d, Sample value %hi"), NumFramesAvailable, i, chunk.chunk[i]);
        }
    }

    if (nonZero)
    {
        m_queue.push(chunk);
    }

    return 0;
}

This doesn't seem to be very good code, and it certainly isn't well commented to explain why it is doing what it is doing.

To answer your question at the "*****  WHAT'S GOING ON HERE  *******" mark:

1) The "* 2" in the chunk.size calculation is clearly because the data is assumed to be stereo.

2) Each sample value is a 2 byte int (16 bit), the chunk needs to have chunk.size samples.

3) The multiplier is to compensate for the memcpy working differently on different platforms.  A char may be one or two bytes,  memcpy moves char sized pieces.  It used to be strictly bytes, but that isn't the case with this code.  This problem would be more efficiently solved with a global DEFINE.

4) The memcpy copies the data from the input buffer to the storage location

5) The loop looks for a silent chunk, allowing a little tolerance -1 to 1 is silent.

Now, here is why I think this is poor code:

For ordinary audio, there are going to be a heck of a lot of Log messages generated.
If it is non-zero (not silent) the chunk will be kept (pushed onto the queue).  If logging every value isn't needed, the loop should break as soon as the chunk is determined not to be silent
The logic of zeroing 1s and -1s in a kept chunk makes no sense to me
If the chunk is not kept, the chunk.chunk array is not deleted (assuming the kept chunks are deleted downstream).  This is a memory leak.

That's enough for now, hope it was helpful.

Ced


SMA, SVM and PSD python I'm looking for python functions (package) to calculate SMA (signal magnitude area), SVM (signal vector magnitude) and PSD (Power Spectral Density). The goal is to extract features from an accelerometer data.
Thanks

You can use matplotlib.pyplot.psd to get PSD (Power Spectral Density): it returns the values for the power spectrum and the corresponding frequencies.
As for SMA (Signal Magnitude Area) and SVM (Signal Vector Magnitude), I am not aware of their support in existing packages. They seem to be easily implementable in one-line functions though.


What is the shortest time period that can be used to find the frequency of a signal burried in noise and harmonics? I have a continuous-time signal that is of the form
$$x(t)=A \sin(2 \pi f t+\phi_1)+B \sin(2 \pi (5f) t+\phi_2)+C \sin(2 \pi (7f) t+\phi_3)+n(t)$$
sampled at 50us where $f$ is the frequency I want to find out A,B and C are constants where $$0.15 \cdot A>B+C+|n(t)|$$. I don't need to know the value of any of the phases or amplitudes. I also know that the value of $f$ is 50Hz +/-15 Hz though i want to know it to at least +/- 1 Hz accuracy.

What ive tried so far:

A PLL-because $f$ is quite a low frequency and my frequency range is so wide and harmonics were present, the bandwidth of my loop filter was quite limited, I couldn't beat a settling time of 100ms(within 1 HZ).
DFT-Sampled over 100ms, the spectra was only really clear when $f$=50Hz.
Measuring the time between zero crossings(Though this falls apart when the amplitudes of the harmonics and noise gets quite high)

Code used for DFT

T = 0.1; % sample time 100ms
Ts = 50e-6; % Sampling period
fs = 1/Ts; % Sampling rate
t = 0:Ts:T-Ts; % Time vector of length 100ms
f1=50;
f2=35;
f3=65;
x1 = sin(f1*2*pi*t)+0.1*sin(3*f1*2*pi*t+pi/3)+0.1*sin(5*f1*2*pi*t+pi/3)+0.01*sin(50*f1*2*pi*t+pi/3);
x2 = sin(f2*2*pi*t)+0.1*sin(3*f2*2*pi*t+pi/3)+0.1*sin(5*f2*2*pi*t+pi/3)+0.01*sin(50*f2*2*pi*t+pi/3);
x3 = sin(f3*2*pi*t)+0.1*sin(3*f3*2*pi*t+pi/3)+0.1*sin(5*f3*2*pi*t+pi/3)+0.01*sin(50*f3*2*pi*t+pi/3);
N = length(t); % Number of elements in the time vector
F = fs/(N); % Frequency step for FFT display
f = (-fs/2):F:(fs/2)-F; % Frequency vector(same length as the time vector)

X=fftshift(fft(x1))/N; % FFT of the signal x1  (blue)
X2=fftshift(fft(x2))/N; % FFT of the signal x2 from task 2 (red)
X3=fftshift(fft(x3))/N; % FFT of the signal x3 from task 3 (black)

figure(2);
subplot(2,1,1), plot(t, x1, 'b', t, x2, 'r', t, x3, 'k'), axis([0 .1 -1.5 1.5]);
legend('x=50Hz','x=35Hz','x=65Hz')
subplot(2,1,2), plot(f, abs(X), 'b', f, abs(X2), 'r', f, abs(X3), 'k'), axis([-100 100 0 0.5]);

if you know that the 5th and 7th harmonic are precise harmonics of the fundamental at frequency $f$, you can use them along with the fundamental to help against the noise.  (the signal is not "buried" in its own harmonics, the harmonics are part of the signal, unlike the noise.)

this is the "pitch detection" problem (even if you're only interested in frequencies around 50 Hz).  perhaps considering computing the Average Squared Difference Function (ASDF)

$$ Q(L) = \frac{1}{M} \sum\limits_{i=0}^{M-1} \left(x[n-i] - x[n-L-i] \right)^2 $$

try computing that for lags $L$ around the period length ($\frac{1}{50}$ second).  pick the $L$ such that $Q(L)$ is minimum.  the noise will keep it from getting to zero, but it will get close to zero.

as to the question in the title of this thread, the shortest time to estimate the frequency is about 1 period ($\frac{1}{50}$ second) plus a little more.


How to write lowpass filter for sampled signal in Python? I have some signal that sampled each 1 ns (1e-9 sec) and have, let say, 1e4 points. I need to filter high frequencies from this signal. Let say I need to filter frequencies higher than 10 MHz. I want that for frequencies lower than cutoff frequency signal will be passed unchanged. It means gain of the filter will be 1 for frequencies lower than cutoff frequency. I would like to be able to specify filter order. I mean, first order filter have 20 db/decade slope (power roll off) after cutoff frequency, second order filter have 40 db/dec slope after cutoff frequency and so on. High performance of code is important.

The frequency response for the filter designed using the butter
function is:

But there is no reason to limit the filter to a constant monotonic
filter design.  If you desire a higher attenuation in the stopband
and steeper transition band, other options exist.  For more information
on specifying a filter using iirdesing see this.  As shown by
the frequency response plots for the butter design the cutoff
frequency (-3dB point) is far from the goal.  This can be alleviated
by down-sampling before filtering (the design
functions will have a difficult time with such a narrow filter, 2% of
the bandwidth).  Lets look at filtering the original sample
rate with the cutoff specified.

import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

from scipy.signal import fir_filter_design as ffd
from scipy.signal import filter_design as ifd

# setup some of the required parameters
Fs = 1e9           # sample-rate defined in the question, down-sampled

# remez (fir) design arguements
Fpass = 10e6       # passband edge
Fstop = 11.1e6     # stopband edge, transition band 100kHz
Wp = Fpass/(Fs)    # pass normalized frequency
Ws = Fstop/(Fs)    # stop normalized frequency

# iirdesign agruements
Wip = (Fpass)/(Fs/2)
Wis = (Fstop+1e6)/(Fs/2)
Rp = 1             # passband ripple
As = 42            # stopband attenuation

# Create a FIR filter, the remez function takes a list of
# "bands" and the amplitude for each band.
taps = 4096
br = ffd.remez(taps, [0, Wp, Ws, .5], [1,0], maxiter=10000)

# The iirdesign takes passband, stopband, passband ripple,
# and stop attenuation.
bc, ac = ifd.iirdesign(Wip, Wis, Rp, As, ftype='ellip')
bb, ab = ifd.iirdesign(Wip, Wis, Rp, As, ftype='cheby2')

As mentioned, because we are trying to filter such a small percent
of the bandwidth the filter will not have a sharp cutoff.  In this
case, lowpass filter, we can reduce the bandwidth to get a better
looking filter.  The python/scipy.signal resample function can
be used to reduce the bandwidth.

Note the resample function will perform filtering to prevent
aliasing.  Prefiltering can also be perfomed (to reduce aliasing)
and in this case we could simply resample by 100 and be done, but
the question asked about creating filters.  For this example
we will downsample by 25 and create a new filter

R = 25;            # how much to down sample by
Fsr = Fs/25.       # down-sampled sample rate
xs = signal.resample(x, len(x)/25.)

If we update the design parameters for the FIR filter the
new response is.

# Down sampled version, create new filter and plot spectrum
R = 25.             # how much to down sample by
Fsr = Fs/R          # down-sampled sample rate
Fstop = 11.1e6      # modified stopband
Wp = Fpass/(Fsr)    # pass normalized frequency
Ws = Fstop/(Fsr)    # stop normalized frequency
taps = 256
br = ffd.remez(taps, [0, Wp, Ws, .5], [1,0], maxiter=10000)

The filter operating on the downsampled data has a better
response.  Another benefit of using a FIR filter is that you
will have linear phase response.


Alternative function for Matlab Stationary Wavelet Transformation I'm currently working on a project of mine where I have to use stationary wavelet transformation. I found a perfect way to do the transformation using wavelet toolbox of the Matlab. But since it is not free, I need to find a free function to do the exactly same thing SWT function of wavelet toolbox does.

What I Found

I found a Matlab functions package called WaveLab. They have defined a function called FWT_Stat to do Stationary wavelet transformation.

FWT_Stat -- stationary wavelet transform
Usage
 StatWT = FWT_Stat(x,D,qmf)
Inputs
 x        array of dyadic length n=2^J
 L        degree of coarsest scale
 qmf      orthonormal quadrature mirror filter
Outputs
 StatWT   stationary wavelet transform table
          formally same data structure as packet table
          log_2(n)-D scales by n elements

Below shows information about SWT function from wavelet toolbox.

SWC = swt(X,N,'wname') computes the stationary wavelet decomposition of the signal X at level N, using 'wname'.

Problem

FWT_Stat function asks for 3 parameters and final parameter is orthonormal quadrature mirror filter.But I couldn't find a way to add haar filter as the input for 3rd parameter .If I used SWT function from wavelet toolbox, all I have to do is add wavelet name (e.g: harr) instead of wname and swt function take care of all other things.

So my problem is, how can I add haar filter as the 3rd parameter for FWT_stat function on Matlab without Wavelet toolbox.

Or, Is there any alternative function (which is free) to do the exact same thing as SWT function from wavelet toolbox?

I found an answer to this problem and hope to share it here.

Luckily there is another Matlab function call "MakeONFilter" which included in WaveLab package. That function can use to make a filter which can use as the input for FWT_Stat function.

E.g:

haarF = MakeONFilter('Haar');
FWT_Stat(x,D,haarF);


Confusion about shannon's entropy of the grayscale image after lossy compression I have a .pgm format 480 * 640 gray scale image named 'columns.pgm'. Using PCA (principal component analysis), by preserving 40 principal components I compressed original image with 1218 KB to 301 KB using a Matlab function which I wrote. Now about the entropy. The entropy of the original file was 6.93 and of the compressed file was 7.17. Actually the entropy of the compressed version must be lesser, right? I am confused. Can someone give me a intuition based explanation. And are the cases similar for lossy and lossless compression?

Shannon's entropy only makes sense when you define it over a probability distribution. You're using a probability distribution where each pixel is sampled independently from the histogram of pixel intensities. Clearly, real pictures don't come from such a distribution. If they did, they would all look like white noise.

Also note that individual elements don't have entropy. Only distributions have entropy. So if you have two images it doesn't make sense to compare their entropy. It only makes sense to compare the entropy of the random distributions that you sampled the images from.


Why are the real- & imaginary part of complex white gaussian noise independent processes? Assume we transmit a bandpass signal over an AWGN channel adding the Gaussian noise contribution $n$

\begin{equation}
Acos(\omega_ct) + n(t)
\end{equation}

Further, the bandpass signal is generated by IQ-modulation. If we perform homodyne IQ demodulation + lowpass filtering at the receiver should the resulting in-phase & quadrature component not be affected by the same lowpass noise contribution?
Instead the equivalent complex lowpass noise is modeled as

\begin{equation}
N_{LP}=X+jY
\end{equation}
where the imaginary & real part are both Gaussian distributed and independent processes. Hence the real & imaginary part of the transmitted symbol are affected by different noise contributions.

Why is that?

It is simply because each sample of $n(t)$ has a random magnitude and phase by definition given as $n(t) = |n(t)|e^{j\phi(t)}$.  With real and imaginary components as follows:

$$|n(t)|e^{j\phi(t)} =|n(t)|\cos(phi(t))+j|n(t)|\sin(phi(t)) $$

Real: $I(t) = |n(t)|\cos(\phi(t))$
Imag: $Q(t) = |n(t)|\sin(\phi(t))$

Since the phase and magnitude are independent for each sample, the I and Q components will also be completely independent. Given any I value, there is no constraint or dependence on what the Q value can be for that sample, and vice versa. (A dependence would exist if given the phase, or given the magnitude. And given both then there would be a one-one mapping between I and Q)

Consider the opposite case if the I and Q components of a Gaussian Noise process were dependent such as I = kQ, the resulting noise would stay on a fixed angle passing through the origin (such as staying on the 45° line if I = Q) rather than adding a random magnitude and phase to each sample.


Relation between CTFT and DTFT I need to transform this function:

$$
x(t) = 4\sin(20\pi t ) - 5\cos(24\pi t ) + 3\sin(120\pi t )
$$
into a sequence $x(n)$ given that the sampling frequency should be 50 Hz. So that means the sampling interval $T = \frac{1}{50} = 0.02\ s\ $. So would the sequence be:

$$
x[n] = 4\sin(20\pi0.02n ) - 5\cos(24\pi 0.02n ) + 3\sin(120\pi 0.02n)
$$

I dont know how to compare the CTFT of $x(t)$ and the DTFT of $x(n)$. Any idea ?

Your sequence simply is a sampled version of the continuous signal $x(t)$. So

$$x_d(n)=x(nT),\quad T=\frac{1}{50}$$

The relation between the CTFT of $x(t)$ and the DTFT of $x_d(n)$ is that the DTFT is the sum of shifted versions of the CTFT. The spectra are shifted by the sampling frequency. That's the reason for aliasing.


Doppler frequency when down-converting a receiving signal Suppose we have a monostatic radar system, with only one target crossing at a speed $v$ at the main lobe of the antenna. The Doppler frequency for this case will be given by

$
f_d = \frac{2v}{\lambda} = \frac{2v \cdot f}{c}
$

With $f$ being the carrier frequency of the received signal. However, if we downconvert, $f = 0$, which would give a $f_d = 0$, which does not sound right to me. So, how does downconversion affect the Doppler signature?

When there is relative motion between the transmitter and receiver the frequency seen by the receiver is not $F_c$, so let's say you incur a Doppler shift of $\Delta_f$, then the shifted signal frqeuency is $F_c - \Delta _f$ now when you downconvert with $F_c$ the frequency offset of $\Delta_f$ still remains in the signal. That means the baseband spectrum "centred around freqeuency" of $0$ will be shifted by this value

In the receivers, this offset is estimated by the baseband signal processing (by performing correlation of time domain samples, linear shift in freqeuency domain means phase shift in time domain, so these correlations of samples simply leave the phase corresponding to the frequency offset due to doppler), there are other methods as well. Once these are estimated the clock at the receiver is further corrected with the estimation to tune the receiver to the "seen" signal freqeuency. So that the baseband spectrum is again centred around $0$ frequency


Trouble with time-syncing two radio signals using cross-correlation I posted a similar question a while ago here (I'm posting this follow-up since I'd like to focus more on cross-correlation now). I have a setup where I have two software-defined radios connected to different antennas, and a circuit switches the antennas on and off at the same time. I am trying to time-sync the two signals as precisely as possible. The radios are sampling at 2.4 MS/s and I'm sampling about 0.2 seconds of data. Even though I've synced the two clocks of the radios, due to slight hardware differences there will always be a slight delay in one of the signals which is why I need to sync them in software. I have a system that is able to sync the signals most of the time just by getting the signal's envelopes and identifying where the antenna switching happens, but this doesn't work so well when the "silence" part is too noisy / too high amplitude, so I would like to use cross correlation.

Here are two of the signals I'm trying to align: s0, s1. I believe the correct delay of s0 should be approximately 1069.

I'm normalizing both signals to have a max amplitude of 1, and I've also tried downsampling to make it lower detail. I'm just doing signal.correlate(s0, s1) and plotting the result, but I don't even see any local maximum at the correct delay (I know the correct delay from the other method).

Here's the plotted correlation:

And here it is zoomed in to where I'd expect the actual delay to be:

There doesn't seem to be any local max around 1069.

Is anyone able to successfully sync these two signals using cross correlation, or is it just not possible for this kind of data? Everyone doing similar project to me uses cross correlation to do the time synchronization, so I'm not sure why it's not working at all for me (even when I use their code). Any help would be greatly appreciated!

It will work when you take the 2nd gradient of the signals:

import numpy as np
from scipy import signal

s0 = np.gradient(np.gradient(s0))
s1 = np.gradient(np.gradient(s1))
np.argmax(signal.correlate(s0, s1)) -> 525358

That corresponds to a shift of 1071 which is close to your expected 1069
Interestingly the minimum (most negative correlation) is close by at 1069

Plotting the cross-correlation of 2nd gradient gives a reasonably sharp peak:

Zomming in to the accordingly shifted data shows that this worked:

Using the most negative correlation results in the almost same picture:

What might be surprising is that the min and max are just two samples apart and also that the signal seems to be significantly phase shifted. On the other hand the start is matched perfectly. My interpretation of this is that the sync is not so much done by what we see as signal but by some high frequency parts of noise received in both channels. Since this 2 sample shift between min and max is the typical result when matching two channels that received the same random noise.

I can not explain easily why the correlation of the original data works in some cases. But here is an simple example that gives an idea:

a = np.array([0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0, 0, 0, 0], dtype=float)
b = np.array([0, 1, 2, 3, 4, 3, 2, 1, 0, 1, 2, 3, 4, 3, 2, 1, 0, 0, 0, 0, 0], dtype=float)
np.argmax(signal.correlate(a, b)) -> 16

which is correct

but when adding an offset it fails

a += 5
b += 5
np.argmax(signal.correlate(a, b)) -> 20

the correlation of the first gradient still gives the right result but it will fail also if a significant ramp is added. Though for the 2nd gradient the correct result is achieved

a += np.array(range(0,len(a)))*5
b += np.array(range(0,len(b)))*5
np.argmax(signal.correlate(a, b)) -> 20
a = np.gradient(a)
b = np.gradient(b)
np.argmax(signal.correlate(a, b)) -> 20
a = np.gradient(a)
b = np.gradient(b)
np.argmax(signal.correlate(a, b)) -> 16


Libraries for Voice Activity Detection (Not Speech Recognition) As follow up to my previous question I was wondering if there are any speech detection libraries in existence. By speech detection I mean passing in an audio buffer and getting back an index of where speech starts and stops. So if I have 10 seconds of audio sampling at 44kHz, I would expect an array of numbers such as:

44000
88000
123000
190334
...

This would indicate for example that speech starts one second in and then finishes at the two second point, etc.

What I'm not looking for is speech recognition which writes out text from spoken word. This unfortunately is what I see a lot of when I google 'speech detection'.

It would be great if the library was in C, C++ or even Objective-C as I'm writing an app for the iPhone.

Thanks!

In my answer to your that question, I had mentioned that Voice Activity Detection is a standard feature for codecs like G.729 and such others.

You should look for reference encoders and decoders for algorithms that applies this.

One such example is - http://www.voiceage.com/openinit_g729.php

Another possible source is Speex codec. Which implements VAD

BTW: You should google "Voice Activity Detection" or "Talk Spurt" rather than "Speech Detection".


SigLib FSk Function: SDA_FskModulate I would like to implement the example code from SigLib's example "fsk.c".
The defines are below

#define SAMPLE_LENGTH                   ((SLArrayIndex_t)512)

#define NUMBER_OF_LOOPS                 4L

#define SAMPLE_RATE                     9600.0

#define BAUD_RATE                       1200.0
/* Number of samples per symbol */

#define SYMBOL_LENGTH                   ((SLArrayIndex_t)(SAMPLE_RATE / BAUD_RATE))

#define MAX_RX_STRING_LENGTH            ((SLArrayIndex_t)80)
/* Frequency of sine wave in table */

#define CARRIER_TABLE_FREQ              ((SLData_t)100.0)
/* Number of samples in each of cos and sine table */

#define CARRIER_SINE_TABLE_SIZE  ((SLArrayIndex_t)(SAMPLE_RATE / CARRIER_TABLE_FREQ))
/* Must be an integer number of cycles */
/* Rx Filter length */

#define RX_FILTER_LENGTH    ((SLArrayIndex_t)((2 * SYMBOL_LENGTH) + 1))
/* Low carrier frequency */

#define FREQ_1300_HZ                    1300.0
/* High carrier frequency */

#define FREQ_2100_HZ                    2100.0
/* Carrier frequency for '0' */

#define CARRIER_FREQ_ZERO               FREQ_2100_HZ
/* Carrier frequency for '1' */

#define CARRIER_FREQ_ONE                FREQ_1300_HZ
/* Bandwidth of detection filter */

#define FILTER_BANDWIDTH                400.0

The initialization call to the function is this:

`SIF_FskModulate (pCarrierTable, (CARRIER_TABLE_FREQ / SAMPLE_RATE), CARRIER_SINE_TABLE_SIZE);`

In the code, the second parameter: (CARRIER_TABLE_FREQ / SAMPLE_RATE) is commented with:
/* Carrier frequency */

This number works out to a fractional amount - 100.0/9600 = small. How is this used? What would the FSK carrier be set to, or is that determined by the rate at which you load the pData array into the DAC?

Can someone clarify what that means in terms of using this code?

Thank you!

The comment is mistaken.  The carrier frequency is 100 Hz.  The sample rate is 9600 Hz.  The carrier frequency divided by the sample rate is the carrier phase increment (phase change per sample), not the carrier frequency.  The units are radians / $2\pi$ (i.e. the phase increment is normalized from 0 to 1, not 0 to $2\pi$).  That phase increment can be used to index into a sine table to create a sinusoid.


Finding the linear prediction coefficients for a sampled sinusoid I am unable to work my way through this problem -

A sinusoid with random phase is given as () = 2 sin(0.25π + ). If this signal is sampled at 10
KHz, estimate the coefficients of a second order LP model. Find the magnitude and phase values of the
second order LP error filter at  = 0.5. Estimate the bias (error) in the peak frequency (frequency at the
maximum magnitude response) of ()

I am thinking of proceeding with Yule-Walker equations but unable to solve this.

Can you please provide a way to solve this.

In the absence of more details, I assume that

the phase $\Theta$ is a random variable and
$\Theta$ is uniformly distributed on the interval $[0,2\pi]$, so that its probability density function is $\frac{1}{2\pi}$ on that interval.

This is a common choice of model for such problems, and I will show that its guarantees that the stochastic process $x$ is wide-sense stationary (WSS).

\begin{equation}
\begin{split}
\mathsf{E}[x[n]] &=~ \mathsf{E}[2\sin(0.25\pi n + \Theta)]\\
&=~ 2\mathsf{E}[\sin(0.25\pi n + \Theta)]\\
&=~ 2\int_{0}^{2\pi}\sin(0.25\pi n + \theta)\frac{1}{2\pi}d\theta\\
&=~ \frac{1}{\pi}\int_{0}^{2\pi}\sin(0.25\pi n + \theta)d\theta
\end{split}
\end{equation}
The integral is the integral of sine over a $2\pi$-wide interval, so it is equal to 0. Hence $\mathsf{E}[x[n]] = 0$, independent of $n$.

\begin{equation}
\begin{split}
\mathsf{Var}[x[n]] &=~ \mathsf{E}\left[\left(x[n] - \mathsf{E}[x[n]]\right)^2\right]\\
&=~ \mathsf{E}\left[\left(x[n] - 0\right)^2\right]\\
&=~ \mathsf{E}\left[x[n]^2\right]\\
&=~ \mathsf{E}\left[4\sin^2(0.25\pi n + \Theta)\right]\\
&=~ 4\int_{0}^{2\pi}\sin^2(0.25\pi n + \theta)\frac{1}{2\pi}d\theta\\
&=~ \frac{2}{\pi}\int_{0}^{2\pi}\sin^2(0.25\pi n + \theta)d\theta\\
&=~ \frac{2}{\pi}\times\pi ~=~ 2,
\end{split}
\end{equation}
which is independent of $n$. We see that $x$ is a WSS process.

This is important because we rely on expected values of the form
\begin{equation}
\begin{split}
\mathsf{E}[x[n]x[m]] &=~ \mathsf{E}[4\sin(0.25\pi n + \Theta)\sin(0.25\pi m + \Theta)]\\
&=~ \frac{4}{2\pi}\int_{0}^{2\pi}\sin(0.25\pi n + \theta)\sin(0.25\pi m + \theta)d\theta
\end{split}
\end{equation}
I leave it as an exercise to show that the value of this integral depends only on the difference $n-m$ and not on $n$ and $m$ separately.

In particular, we can define $R[i] = \mathsf{E}[x[n]x[n-i]]$ and know that $R[i]$ is really a function of $i$ alone and has no dependence on $n$.

For second-order linear predictive coding (LPC), we estimate $x[n]$ from the previous two values:
\begin{equation}
\widehat{x}[n] = a_1x[n-1] + a_2x[n-2].
\end{equation}
The error is then
\begin{equation}
\begin{split}
e[n] &=~ x[n] - \widehat{x}[n]\\
 &=~ x[n] - a_1x[n-1] - a_2x[n-2].
\end{split}
\end{equation}
The coefficients $a_1$ and $a_2$ must be chosen to minimize mean square error, which is $\mathsf{E}[e[n]^2]$.

To minimize mean square error, the error must be stochastically orthogonal to the random variables that make up the estimate:

\begin{eqnarray}
\mathsf{E}[e[n]x[n-1]] &=& 0,\\
\mathsf{E}[e[n]x[n-2]] &=& 0.
\end{eqnarray}
These will give us linear equations whose solution reveals the proper choice of $a_1$ and $a_2$.

\begin{equation}
\begin{split}
\mathsf{E}[e[n]x[n-1]] &=~ \mathsf{E}\left[\left(x[n] - a_1x[n-1] - a_2x[n-2]\right)x[n-1]\right)]\\
&=~ \mathsf{E}[x[n]x[n-1]] - a_1\mathsf{E}[x[n-1]^2] - a_2\mathsf{E}[x[n-2]x[n-1]]\\
&=~ R[1] - a_1R[0] - a_2R[1]
\end{split}
\end{equation}
Since we want this quantity to be equal to zero, we have the linear equation
\begin{equation}
R[0]a_1 + R[1]a_2 = R[1].
\end{equation}

\begin{equation}
\begin{split}
\mathsf{E}[e[n]x[n-2]] &=~ \mathsf{E}\left[\left(x[n] - a_1x[n-1] - a_2x[n-2]\right)x[n-2]\right)]\\
&=~ \mathsf{E}[x[n]x[n-2]] - a_1\mathsf{E}[x[n-1]x[n-2]] - a_2\mathsf{E}[x[n-2]^2]\\
&=~ R[2] - a_1R[1] - a_2R[0]
\end{split}
\end{equation}
Since we also want this quantity to be zero, we have the linear equation
\begin{equation}
R[1]a_1 + R[0]a_2 = R[2].
\end{equation}

We now have the system of equations
\begin{equation}
\left(\begin{array}{cc}
R[0] & R[1]\\
R[1] & R[0]
\end{array}\right)
\left(\begin{array}{c}a_1\\a_2\end{array}\right)
=
\left(\begin{array}{c}R[1]\\R[2]\end{array}\right),
\end{equation}
whose solution is
\begin{equation}
\left(\begin{array}{c}a_1\\a_2\end{array}\right)
= \frac{1}{R[0]^2 - R[1]^2}
\left(\begin{array}{c}
R[0] & -R[1]\\-R[1] & R[0]
\end{array}\right)
\left(\begin{array}{c}R[1]\\R[2]\end{array}\right)
\end{equation}

I leave it as an exercise to compute the correlations $R[0]$, $R[1]$, and $R[2]$.

An aside: In higher-order LPC we take advantange of the very nice structure of the matrix in the linear equations (it's a Toeplitz matrix). That structure allows one to use the so-called Levinson-Durbin algorithm to solve it numerically. While the matrix is real symmetric and thus is numerically amenable to numerous algorithms, Levinson-Durbin is the one used in practice because of its speed and because it yields other useful quantities as it runs. These reflection coefficients (RCs) that it yields along the way, are proxies for the linear prediction coefficients (LPCs).

I should note here that there are sometimes differences in sign. See, for example, the help page for MATLAB's lpc function. But once the signs are chosen, there is consistency afterward. Check the definition of LPC in your example before starting.

The filter to compute the error is $x[n] - a_1x[n-1] - a_2x[n-2]$, so the impulse response is $\mathbf{h} = (h_0, h_1, h_2) = (1, -a_1, -a_2)$. The corresponding frequency response is
\begin{equation}
H(e^{i\omega}) = 1 - a_1e^{-i\omega} - a_2e^{-2i\omega}.
\end{equation}

In some applications (such as vocoders), the error ($x[n] - \widehat{x}[n]$) is assumed to be a noise process. If the LPCs have been received, then a realization of the noise process is generated, and an approximation of the original signal is created by using the inverse filter, which is an IIR filter with $z$-transform
\begin{equation}
\frac{1}{H(z)} = \frac{1}{1 - a_1z^{-1} - \cdots - a_kz^{-k}}.
\end{equation}


What is the mathematical interpretation of using direct inversion around a single pole/zero Let $F[z]=N^2\frac{z(z-(3+j\sqrt{7})/2)(z-(3-j\sqrt{7})/2)}{(z-1)^3}$

This has 3 poles at $z=1$; one zero at $z=0$; and a conjugate pair of zeros at $z=\frac{3\pm{}j\sqrt{7}}{2}$

Assuming a contour $H$ encircles these roots and zeros, a unique inverse is obtained from direct inversion:

$$f=\frac{1}{2{\pi}j}\oint_H Fz^{n-1} \,dz$$

What, if any, is the mathematical interpretation of taking a contour around a single pole or zero?

For the integral to be an inverse Z transform, the contour $H$ must enclose the origin, and it must be inside the region of convergence (ROC). In this way, the contour defines the ROC. E.g., if it encircles all poles, the corresponding sequence is causal. If all poles are outside the contour, then the resulting sequence is anti-causal. Finally, if there are poles inside and outside the contour, the resulting sequence is two-sided.

In your example, all poles are at $z=1$, so there are only two options: either the contour encloses these poles, then the sequence is causal. Otherwise, the sequence is anti-causal. Note that the location of the zeros is irrelevant here.

If you had two distinct poles, you would have three choices:

both poles are outside $H$: the resulting sequence is anti-causal (left-sided)
$H$ encloses only 1 pole: the resulting sequence is two-sided (non-causal)
$H$ encloses both poles: the resulting sequence is causal (right-sided)


Advantages and disadvantages of multiplying sparse vector with sparse matrix Assume I have a sparse vector $x$ whose few values are non-zeros, that vector is multiplied with sparse unitary matrix $Y$ giving $z = Yx$. Are there any advantages of detecting the sparse vector $x$  based on $z$ the in that case?

You are referring here to sparse signal sampling and Reconstruction. You do not require the sampling matrix to be sparse. Infact we require the sampling matrix to follow the restrictive isometry property. And sparse matrices do not follow this property. Random matrices do. So no point having a sparse sampling matrix


How would I go about detecting and following edges in time lapse data? I have a time lapse movie from a specimen and would like to trace the cell wall movements, but I'm not sure how to. By eye, following the cell wall is fairly easy but it seems computationally quite difficult. I was looking into Caffe and deep learning networks, but the tutorials seem more about finding a pattern in images, not so much about tracing data within a movie. If someone can provide me some pointers on how to go about my data that would be great.

Here is a link to an example movie: Example Movie

I'm only interested in the movements of the outer wall, as indicated here:

Deep learning might have something to offer in the long term but at the expense of training the network with very detailed information about what constitutes an "outer wall". A deep learning network could then iterate all frames and markup those blocks that are very likely to constitute the "outer wall".

In the short term, a technique based on normalised cross correlation might offer more satisfactory results. The idea in this case is simple: Define an area of the image (let's call it a Region of Interest (ROI)) that you want to track and then try to "find" that area within the next frame and the next frame and so on.

The primary assumption here, of course, is that the ROI doesn't change significantly from frame to frame. And in this particular case, this takes us back to what exactly does an "outer wall" mean?

Here is an example of tracking a point on the "outer wall":

And here is the data of the track (frame num, X, Y).

This is the best tracking I could do out of approximately 8 attempts (I will come back to this) and it was done using Blender.

Blender is using the normalised cross correlation method to track "markers" from a video feed. Depending on the type and number of markers tracked, you can do a number of things such as stabilise video or solve for the camera movement and position in order to synchronise a "virtual" and a real scene.

In this particular case, I am only tracking Location. I could track Location, Rotation, Scale or any combination of those but I figured that this phenomenon is happening more or less on a narrow depth of field. So, just tracking for location will not distort the results.

The tracked coordinates are the green and red timeseries, green for a "Y" coordinate and red for a "X" coordinate.

As you can see, it's not a simple linear movement. The tracked marker seems to be prescribing a tiny little ellipse around its position.

And this brings us to the main "problem" you are going to have to deal with: WHAT is an outer wall? Here are a few observations:

The "outer wall" displacement varies across it's length. The "outer wall" of this membrane (?) is roughly parallel to the width of the video. So, the recorded displacement, really depends on WHERE you are going to track. If I tried to track somewhere else along the length of the track I would have had a different profile.
The "WHERE are you going track?" question is made harder, because the "outer wall" is semi-transparent and furthermore, its transparency throughout it's "pumping" cycle seems to depend on the phase of it. So, there are points in the cycle that there are very clearly defined edges (which is good for accurate tracking) but gradually, these points vanish and introduce drift to the tracking. Gradually, the marker will "slip" from the target it is supposed to be tracking and start tracking a slightly adjacent image patch that happens to look similar.

Therefore, what I would suggest is to set-up a similar script that opens a video, establishes a marker patch and then tracks the patch frame-by-frame, using information that is readily available in many different scientific platforms (for example, for another example).

BUT!, having now established tracking, you will have to go ahead and manipulate the successive frames of the video in order to preserve the features. This could be a simple contrast enhancement with a steep window, or a histogram equalisation or a normalisation step or possibly the addition of a high pass filter to make the edges more pronounced and improve tracking.

Finally, if you would like to replicate this tracking, using your video, you will need:

Blender
Some demo on tracking. This one is very good.
Specific parameters

This tracking was produced with: Tracking Settings(Preset, Blurry Footage, Motion Model: Loc, Match: Previous Frame). You can match either on a frame-by-frame mode, or to a keyframe, in which case, you get displacements relative to that particular frame. Might be useful for you if you want to measure everything with respect to a particular point in the cycle.

A script to extract the tracking data to a CSV file. Like this one.

(Of course, you can track as many markers as you like and their data will be put in a separate file).

Hope this helps.


Are Total Harmonic Distortion and Noise (THD+N) and THD correlated in practice? I've read about THD and THD+N due to interest in audio quality.

Some audio devices (ej. speakers) report %THD or %THD+N. %THD+N is greater than %THD. I've found that sometimes the difference is huge, which makes me think that a significant part of the noise is due to large generation of spurious frequencies (which are not harmonics of the reference input frequency).

I think that in similar circumstances non-harmonic frequencies should generate to the hearer a perceived distortion greater than those due to harmonics frequencies (Am I right?). If this is true, the %THD+N must be much more relevant than %THD. In such case %THD can only be relevant for comparing two audio devices if it is positively correlated with %THD+N.

Q1 Are %THD+N and %THD correlated in practice in the case of common audio devices (speakers, headphones, amplifiers)?

Q2 How important is %THD in the reproduction of the input signal in speakers or headphones?

Q1 Are %THD+N and %THD correlated in practice in the case of common
  audio devices (speakers, headphones, amplifiers)?

This really depends. THD is just a contributor to THD+N. If the noise dominates the numbers can be vastly different and if the THD dominates the numbers are basically the same. Passive components (passive speakers or headphones) have no inherent electrical noise so THD tends to dominate unless you hit a point where you trigger mechanical or acoustic noise as well. A good example of non-harmonic noise is port noise on a vented speaker.

  Q2 How important is %THD in the reproduction of the input signal in
  speakers or headphones?

Depends on your definition of important on what you care about, but specifically in speaker design we keep a very close eye on THD since speakers are fairly non-linear.


Inverse of the Gibbs Phenomenon to generate a sine wave I am a quite new to the area of digital signal processing. I am currently doing a project on sine wave generation. I know that according to Gibbs phenomenon, many sine waves add up to form a square wave. So, square wave would contain many frequencies. If we filter that square wave stage by stage using low pass and band pass filters, will I be able to get a sine wave at expected fundamental frequency? Please help.

Square waves are composed of sine waves at the odd harmonics- i.e. 1st, 3rd, 5th, etc.  So yes, you can produce a sine wave from a square wave if you filter it with a low-pass filter whose bandwidth is at least as big as the first harmonic frequency, and whose cutoff frequency is lower than the third harmonic.


LTI system impulse response output If I have a LTI system which has impulse response $h(n) = \left(\frac{1}{2}\right)^nu(n)$   and the entry is $x[n]=u[n]-u[n-1]$ where $u[n]$ is the unit step function. In order to find the output I have to convolute the $h[n]$ with $u[n]$ and knowing $x(n)=u[n]-u[n-1]\Longrightarrow y(n)=h(n)+h(n-1)$ take the result. Is my thought right? Because I saw somewhere some examples solving differently.
for the above example the result is $y[n]=\left(\frac{1}{2}\right)^nu(n)-\left(\frac{1}{2}\right)^{n-1}u(n-1)  $
(there is a similar question link but I have confused from the discussion.)

If your input signal is

$$x[n]=u[n]-u[n-1]$$

then it is in fact a unit impulse: $x[n]=\delta[n]$, which has only one non-zero value for $n=0$ ($x[0]=1$). All other values are zero. So the response to this input signal is the system's impulse response, because that's the definition of the impulse response: $y[n]=h[n]$

Note that $u[n]$ is the unit step, and the response to a unit step is not the impulse response but the step response. But for this example you don't need to compute the step response, because, as already mentioned, the input signal is simply a unit impulse.


2-D Distance for Kalman measurements I have designed a simple 3-state $[x\ \dot{x}\ \ddot{x}]$ Kalman filter which is updated by measurements for $x$ and $\dot{x}$. The filter is tracking a peak on a surface in the $x/\dot{x}$ plane. At each iteration, I detect a number of peaks on the surface, each with an $x/\dot{x}$ measurement pair.

I need to select the peak which is closest to the current state as input for the next filter iteration. At the moment I am simply finding the difference in the $x$ and $\dot{x}$ dimensions, scaling them so they are similar in scale, and then computing the 2-D distance. This has been working fine for me, but I know it is not the "right" way to calculate the distance.

Is there a way to use the Kalman variables to determine which measurement is closest to the current estimate?

Not sure where this belongs, so I have this posted in stack, dsp, math, and stats exchanges.

I thought about it a bit more (power of posting), and decided the best way to scale $x$ and $\dot{x}$ is by the standard deviation, which is being tracked within the covariance matrix.

errX = (x_n - x_n-1)/sqrt(var(x));
errXdot = (xdot_n - xdot_n-1)/sqrt(var(xdot));
err2D = sqrt(errX^2 + errXdot^2);


How to scale the amplitude of an FFT of a chirp signal I'm analysing some accelerometer data. The accelerometer was placed on a mechanical shaker and on this shaker a frequency sweep was performed.

When I look at the amplitudes of single sine waves they see to fit in the time and frequency domain. However, if I perform a sweep, the amplitudes are wrong.

The accelerometer data should be scaled right, I think it has to be some MATLAB fft specific stuff. The following code was used to perform the FFT:

nfft= length(signalToAnalyse);
%generate the freqvec
halfn = floor(nfft / 2)+1;
deltaf = 1 / ( nfft / Fs);
freq = (0:(halfn-1)) * deltaf;
% perform FFT
fftRawResult = fft(signalToAnalyse,nfft);
% convert from 2 sided spectrum to 1 sided
magnitude(1) = abs(fftRawResult (1)) ./ (nfft);
magnitude(2:(halfn-1)) = abs(fftRawResult (2:(halfn-1))) ./ (nfft / 2);
magnitude(halfn) = abs(fftRawResult (halfn)) ./ (nfft);

Anybode some Ideas?

Below is the signal in time domain and frequency domain.

First, thanks for all your inputs. I did some research and I found my mistake. The spectrum of a linear chirp can, according to Wikipedia (search for: "chirp spectrum wiki") be analytical described in the frequency domain. The formula for the amplitude

gave me the hint. I just used the constant factor and a factor of samplingfrequency/numberofsamples to match the FFT. I do not know exactly why, but it seems to fit.

Following Picture shows the result:

Following code can be used to simulate it:

    %% Test script
close all
clear all
clc
% create timevector
L =1000 ;
aux = 0:1:(L)-1;  % indexvector
% get fs
fs = 1000;
% create timevector
t = (1/fs)*aux;
f0 = 1;
f1 = 100;
% chirp in timedomain
SignalChirp = chirp(t,f0,max(t),f1);
% Window length
nfft= length(SignalChirp);
%generate the vector of frequencies
halfn = floor(nfft / 2)+1;
deltaf = 1 / ( nfft / fs);
ffft = (0:(halfn-1)) * deltaf;
% perform FFT
X = fft(SignalChirp,nfft);
magni1(1) = abs(X(1)) ./ (nfft);
magni1(2:(halfn-1)) = abs(X(2:(halfn-1))) ./ (nfft / 2);
magni1(halfn) = abs(X(halfn)) ./ (nfft);

figure(1)
subplot(2,1,1)
plot(t,SignalChirp,'k');
title(['Chirp Signal']);
xlabel('Time(s)');
ylabel('Amplitude');
subplot(2,1,2)
plot(ffft,magni1),hold on
duration = max(t);
c = (fs/L);
% https://en.wikipedia.org/wiki/Chirp_spectrum
scalingLinChirp = sqrt(duration/((f1-f0)));
plot(ffft,ones([length(ffft) 1])*scalingLinChirp*c), hold on
xlabel('Freq(Hz)');
legend('FFT Result' ,'Calculated Amplitude' )
ylabel('Amplitude');

If I use this scale factor to my FFT of my accelerometer data it matches the calculated theory much better than before. Sadly I can’t show you guys another picture due my low rep.

Thanks for all your Inputs!


Continuous speech recognition with HMM I don't understand how this works. Help would be really appreciated.

Let's say we have a 10 second sound input, we make a feature vector every 10ms - so we have 1000 vectors - $\mathbf{o}=o_1, \cdots,o_{1000}$.

For simplicity, let's say that all we want to know, is what sequence of phones $p=p_1,\cdots,p_n$ is the answer to $\text{argmax } P(o|p)$ for $p$, or at least what's a good estimate/close answer. For all phones, we have a seperate left to right HMM model (I think that's what makes sense in this scenario).

How exactly is this done?
Will similar enough vectors be "put together" and only treated as 1?
The input is probably not one phone or a thousand phones - but we don't know that... so for what $p$ will $P(o|p)$ be computed?

Will similar enough vectors be "put together" and only treated as 1?

No, this is not how it is done.

  so for what p will P(o|p) be computed?

Your assumption that recognition is performed by considering one by one all possible phone or word sequences and scoring them is wrong. For just a few seconds of speech the number of possible words would be extremely large! This approach - scoring each candidate model among the entire space and picking up the one with the largest score, is actually used only for some small vocabulary / single-word applications (say recognizing several options on a voice menu). I think this confusion stems from the fact that the classic Rabiner tutorial on HMM makes a lot of sense for finite-vocabulary applications, but stays away from all the difficulties/technicalities of continuous applications.

In continuous speech recognition there are several layers on top of the phones models. Phone models are joined together to form words models, and words models are joined together to form a language model. The recognition problem is formulated as searching the path of least cost through this graph.

One could do the same by piecing together HMMs, though... Let's take a simpler example... You have 5 phones and you want to recognize any sequence of them (no words models, no language models - let's consider meaningless speech for this example). You train one left-right HMM per phone, say 3 or 4 states - maybe on clean individual recordings of each phone. To perform continuous recognition, you build a big HMM by adding transition states between the last state of each phone and the first state of each phone. Performing recognition with this new model is done by finding the most likely sequence of states given the sequence of observation vectors - with the Viterbi algorithm. You will know exactly which sequence of states (and thus phones) is the most compatible with the observation. Practically, things are never done this way, and this is why I downplayed in my previous reply the importance of Rabiner's "problem 2": in a true continuous application, the size of the vocabulary is such that considering this giant HMM made of all possibles connections between phones to make words, and words to make sentences would make the naive use of the Viterbi algorithm impossible. We stop reasoning in terms of HMM and we start thinking in terms of FST.


Significance of signal energy in audio? Probably a dumb question but I can't find the answer anywhere. It seems to me that the energy carried by a sound wave should be proportional to frequency squared. For example, the waves $x(t)=\sin(t)$, $x(t)=\sin(100\cdot t)$, and $x(t)=\frac{1}{\sqrt{2}}$ carry different amounts of energy - the last one carries no energy at all and sounds like silence. But in DSP they are all considered to have the same signal energy (integral of $x(t)$ squared). What's the significance of that?

It seems to me that the energy carried by a sound wave should be
  proportional to frequency squared.

Why? What makes you think that this is the case ?

  But in DSP they are all considered to have the same signal energy
  (integral of x(t) squared)

The concept of associating power or energy with a a digital signal is somewhat tricky. Digital signals are just a bunch of numbers so assigning actual physical properties to them doesn't really work unless you carefully define the context.

Mostly the concept is used assuming that the digital signal is an accurate representation of an analog signal that has actual physical properties. If this relationship is well defined and constant you can compare the sum-of-squares of two different digital signals and draw some conclusions on the energy difference of the actual analog situations. However, that really requires "all things being equal" and there are actually a lot of things that need to be equal for this to hold.

STEP 1: basic physics of power

Power, intensity or energy are typically defined or calculated as the product of two field quantities (not one). In electricity power is defined as the product of of voltage and current. Simple example: if your signals were defined as "voltage over a one Ohm resistor", than you can indeed calculate the power by summing the squares and it would be the same power for all of your three examples since the current is proportional to the voltage.

If it were the voltage over an ideal inductor or capacitor then the answer would be wrong. Average power over one of these would be zero in all cases. For an ideal inductor voltage and current are 90 degrees out of phase and the average power is therefore zero.

STEP 2: physics of sound

The two field quantities that make up a sound wave are sound pressure and particle velocity. First you need to define which one your digital signal actually represents. Most of the time it's pressure but many second order microphones (cardioid, dipoles) will create a signal that's proportional to the particle velocity. If you are far enough away from the sound source than pressure and velocity are in-phase and proportional to each other and the proportionality factor is the free field impedance of air (density times speed of sound). In this case you can actually calculate the intensity at the microphone position by simply squaring and scaling. If you are close to a sound source that doesn't work anymore and you actually need to measure  both signals.

As others have said: per it's physical definition sound can't exist at 0 Hz. Sound pressure is a variation of pressure around a steady-state average and if there is no variation there is no sound.

Step 3: human perception

Sound energy and perceived loudness are only loosely related. Human perception is very dependent on frequency and its quite non-linear. Building good loudness models is quite difficult. Simple example: 1 Watt of sound at 1kHz radiated into a typical residential room would extremely loud and quite painful. At 50 kHz it would be utterly inaudible, even though it's the exact same amount of physical sound power.


plotting phase of a signal adding delay I'm trying to plot the phase of this signal $s(f)=A^2T^2sinc^2(Tf)e^{-(j\pi Tf)}$
How can I plot manually this signal?I have to follow some particular rules?I have problems with the delay.

Edit

After calculation I found that the phase for this signal is: $\phi(s(f))=-\pi Tf$

it should be this the plot. Right or not?

I plot it for $T =3$

This is indeed correct. There is a few things you could do to make it a better graph:

Label the X-axis
Label the Y-axis
Use real physical units if it's applicable. Phase is measured either in radians or in degrees. These are VERY different things, so proper labeling with units helps to clarify what you are using. Similar if f is a frequency in Hz (or kHz or MHz) and T is a time (in nanoseconds or fortnights), it should be stated.
You could consider changing the Y-grid so that the grid lines match with meaningful phase values such as 2*pi (or an integer multiple thereof).
Along the same lines, you could alternatively normalize the phase to a meaningful value (like 2*pi) and keep the integer grid.
If the result was obtained by Fourier transform of a real valued time domain signal and the context is clear, consider omitting the left half (f<0) of the plot. The convention is to only plot positive frequencies and the complex conjugate symmetry is implied.

This may sound nit-picky, but making a good graph is typically time well spend. It really helps to clarify what's going on, simplifies communication and will help you to remember what you did when you look at it in 3 months from now.


Can somebody help me with some very basic questions about the fundamentals of information theory? I just started to study information content of the source and there are some very fundamental topics I am not really clear on. These might be a long question or a series of questions. Please help me.

Relationship between uncertainty and information:
Information is related to uncertainty or surprise. What is the nature of this relationship? Is there any difference between uncertainty and information? Before the information reaches the destination, there is an uncertainty at the destination about which signal would appear. When the information reaches the destination, there is no uncertainty left, right? So, information is literally killing or eliminating the uncertainty. Is my concept about information correct?
Information content of an improbable or unlikely message:
The relationship between the probability of a message and information content is $I=\log_2(1/P)$, right? So If the probability of occurrence of a symbol is one, then the destination is always sure that that definite symbol will occur. So there is no uncertainty about the message. In other way, when the message reaches the server, there is no surprise. So the information content of that symbol is zero. Mathematically it is zero too. That is, the source didn't really need to transmit that message as the destination knew what message was coming beforehand, correct?

But what would be the information content of a message with probability zero? I know mathematically it is infinity. But what would be the physical interpretation? And the probability of the message to occur is zero, meaning that symbol is never going to be generated. So the destination is never waiting for that symbol to occur. So in that case, isn't the information content of that message supposed to be zero?

That is going to be my first question. I might need to post a few questions later. Thanks in advance.

What is the nature of this relationship? There are many possibilities. Shannon used a logarithm because he wanted information to be additive, that is, $I(xy)=I(x)+I(y)$ for independent events $x$ and $y$.

You're correct that the arrival of a message at the receiver decreases its uncertainty. It may not be decreased to zero, because there may be an error in the communication. Also, a given channel can only support a certain information rate, which is its capacity.
It does not make sense to ask about messages of probability zero. The source alphabet is defined so that all symbols have probability greater than zero.


Converting from stereo to mono by averaging Hi I'm doing a past exam question about converting from stereo to mono.  The scenario is below:

The first question is:  Explain why would this result in a noisy encoding.  How is the noise introduced and what does it sound like?

Second Question:  What would you do to prevent such noise being introduced?

Answer 1:

The only thing I could think of on this is that the two channels could be slightly out of phase and therefore interfere destructively when summed.  I'm not sure what this would sound like but I guess the amplitude and frequency are being changed by this so I know it wouldn't sound right!

Answer 2:

We would have to do something to remove the phase difference, right?

There are indeed rare instances of recordings in which the L&R channels are in opposite phase. This was done in the late 70s and early 80s to create a cheap "spatial sound" effect (like the "surround" button on low-end hi-fi systems of the time), and on such recordings, computing the average of the left and right channels gives a null output. Some stereo effects (like chorus or reverb) would also make a recording sound less lively or flatter after a mono conversion, but this doesn't really qualify as noise. All these are rare problems, and averaging the L&R channels is the standard way of converting from stereo to mono. In the past I've worked on a system that handled large catalogues of audio files, and my approach was to check whether the energy of the average of the channels was below a certain fraction of the energy of the individual channels ; and if so, used the difference instead of the sum as a mono conversion. Ideally you could compute the cross-correlation between the left and right channels to find the lag at which they are maximally correlated, and delay one channel by the lag before the mixing so they are maximally in-phase with each other, but this is overkill for just handling the 15 or so tracks out of 100,000 that had a massive phase inversion problem - to give you an idea of the rarity of the problem...

The most obvious error is converting from 44kHz to 22kHz by skipping every second sample. This will alias all frequencies in the 11kHz - 22kHz range, and will make the downsampled recording sound "metallic", "crunchy" or "gritty". A proper way of downsampling a recording is to first apply a low-pass filter to remove all frequencies above half the target sample rate ; then decimate (skip every second sample).


Translating Matlab scrambler function to C I am trying to implement the following matlab function in c:
Mathworks scrambler Function Documentation

This is the code I've written:

#define N 5
#define MAX 70
#define BITS 7
//char* msg = "This message is Coming to You!";
//char* msg = "Hello world 000\nHello world 001\nHello world 002\nHello world 003\n";
//char* msg = "Hello world 000\nHello world 001\n";
char* msg = "Hel";

uint8_t current[N-1] = {0};
uint8_t to_Xor[N-1] = {0};

uint8_t bin[MAX*BITS];

int main()
{
    int len=strlen(msg);
    int kk = 0;
    int m = 0;
    uint8_t letter = 0;
    uint8_t chars = 0;

    for(int nn = 0; nn < len*BITS; nn++)
    {
        bin[nn] = 0;
    }

    for(int ii = 0; ii < len; ii++)
    {
        letter = msg[ii];
        for(int jj = 6; jj>=0; jj--)
        {
            bin[kk]=(letter&(1<<jj))>>jj;
            //printf("%d\n",bin[k]);
            kk++;
        }

    }

    uint8_t xor = 0;
    uint8_t polynom = 0x0;
    uint8_t scrmbld_output[len*BITS];

    for(int c = 0; c < len*BITS; c++)
    {
        scrmbld_output[c] = 0;
    }

    for(int bin_ndx = 0; bin_ndx < len*BITS; bin_ndx++)
    {

        for(int n = 0; n < N-1; n++)
        {
            if((1<<n)&polynom)
            {
                xor ^= current[n];
            }
        }

        for (int k = N-2; k > 0; k--)
        {
            current[k]=current[k-1];
        }

        current[0] = xor ^ bin[bin_ndx] ^ 1;

        scrmbld_output[bin_ndx] = current[0];
        xor = 0;

        printf("original: %d scrmbld: %d\n",bin[bin_ndx], scrmbld_output[bin_ndx]);
    }

    int k = 0;
    m = 6;

    for(int i = 0; i < len; i++)
    {
        for(int j = 6; j>=0; j--)
        {
            chars |= bin[k]<<m;
            k++;
            m--;
        }
        printf("%c", chars);
        chars = 0;
        m = 6;
    }

    printf("\nlen: %d", len);
    return 0;
}

The results I get from matlab are different.

If I create the following vector and scrambler in matlab:

a = [ 1     0     0     1     0     0     0     1     1     0 ];
scrmblr = comm.Scrambler(2,[1 1 1 0 1], [0 0 0 0]);
ans = scrmblr(a');
ans' = [1     1     0     0     1     0     1     0     1     1]

but what I'm getting is:

original: 1 scrmbld: 0
original: 0 scrmbld: 1
original: 0 scrmbld: 0
original: 1 scrmbld: 0
original: 0 scrmbld: 0
original: 0 scrmbld: 0
original: 0 scrmbld: 1
original: 1 scrmbld: 1
original: 1 scrmbld: 1
original: 0 scrmbld: 1
original: 0 scrmbld: 0
original: 1 scrmbld: 0
original: 0 scrmbld: 1
original: 1 scrmbld: 0
original: 1 scrmbld: 0
original: 1 scrmbld: 1
original: 0 scrmbld: 1
original: 1 scrmbld: 1
original: 1 scrmbld: 0
original: 0 scrmbld: 1
original: 0 scrmbld: 0
original: 0 scrmbld: 0
original: 0 scrmbld: 0
original: 0 scrmbld: 0

I don't see where I went wrong. Any help greatly appreciated, thank you!

I think the mistake would be this line

current[0] = xor ^ bin[bin_ndx] ^ 1;

Remove the xor with 1 at the end.

current[0] = xor ^ bin[bin_ndx];

The MATLAB figure shows 1 because that is the multiplier (scale factor) and may not mean xor with 1.


How to pad a array representing music to pad the stft of this array? I have an array named samples of size (2, 1321967) representing an audio file. When I pass it to scipy stft with

f, t, Zxx = scipy.signal.stft(samples, sample_rate, nperseg=254)

My array t is of size (10411). I would like to pad t and Zxx 127 points to the right and to the left. To get my padding, I would like to pad samples and then pass it to stft.

How many points do I need to pad to samples in order to get the right size for t?

After some reading, I found the solution to my problem. In the default mode, windows overlap on half their size. Scipy starts by adding zeroes at the boundaries on the left and on the right so that the first and last window are centered on the first and last point of data. Then it pads the end with zeroes so that the last windows on the right is full of point. So, in order to add a window on the left/right, you have to add nperseg/2(127 in my case) points of data. So to add 127 windows on the right and on the left, you just need to add 127^2 = 16129 points on the right and of the left of the original data.

Edit : Based on how the stft is computed, if the original signal is zero-padded on the right and left, any new window created will be on data points with value 0, and as a result will have a value of 0. This means that it is exactly equivalent to 0-pad the original signal of 127*n points or to 0-pad the stft of n points.


Sine Wave "Lobes" (Basics, and probably a dumb question) I am beginning my journey on trying to understand a little deeper on signal processing. I wonder if the effect I am seeing is worth any spit or is there a deeper level of understanding is apparent I do not grasp.
I was initially using Jupyter Notebook to make some examples of Aliasing so that I could see the frequency foldback, and give myself a workspace to intuitively understand this concept. And after playing around with drawing sine waves, If I didn't create a sine wave at frequencies based on power 2 of my "sampling rate" then I would create what looked like lobes in my sine wave.
My initial code looks like this
# Aliasing Example
import numpy as np
import matplotlib.pyplot as plot
import math

%matplotlib qt
plot.ion()

# Number of samples
N = pow(2,14) # 16,384 Samples

# Sampling Frequency
Fs = pow(2,14) # 16,384 Samples/second

# Sampling Frequency
T = 1/Fs # 1 second of data

# The amount of sampling time
time = np.linspace(0.0,(N-1)*T,N)

#Nyquist Frequency
Fn = (Fs/2)

# Number of lobes
num_lobes = 4

#Signal 1
f1 = 4096 - num_lobes/4
str1 = "Sine Wave @" + str(f1) + "Hz"
w1 = 2*math.pi*f1
sig1 = np.sin(w1*time)

fig1 = plot.figure()
ax = fig1.add_subplot()
ax.plot(time,sig1)
ax.set_title(str1)
ax.set_xlabel('Time')
ax.set_ylabel('Amplitude')

The plot of the above graph looks something like
It looks like there is an accumulation of phase difference based between my "sampling frequency" of 16384 Hz and the signal frequency 4097 Hz. I don't think I have explored if this actually means anything? Does this have any real-world implications on how I capture data? On a relationship I need to try to satisfy to correctly recreate/capture data? That if I run an FFT or get any other signal information, am I still getting the correct value on the power of the spectrum?
I guess I am looking at this thinking if I were to capture this on a scope I would think I am getting multiple signals going on at once, but it's only a single wave. And that I don't have proper understanding of setting sampling for my signals of concern?

You're just seeing the amplitude of the individual samples, not the wave that travels between them.    If you generate a high-frequency sine wave, you will see that the samples don't necessarily go anywhere close to the actual peaks of the wave.
At some points, it looks like you might expect, with samples near the peaks of the waveform:

but at others, the samples are not near the actual peak of the waveform:

(These are both from a 12.001 kHz wave viewed in Ocenaudio sampled at 48 kHz)
If you are displaying only the samples (as you are in matplotlib because it does only linear interpolation between them), then you will see this as "lobes" that go up and down as the samples move closer and farther from the actual peaks.  But the actual waveform the samples represent does not go up and down, it's constant.


How to convert single output of microphone into multiple values for using fft? i am trying to make a real time audio spectrum display.
My problem is microphone gives an single analog output, but to apply fft we need atleast 8 inputs to fft. how can i convert a single input from microphone to  atleast 8 values for using fft??????

The FFT will process the samples bs time of your analog signal. For 8 consecutive samples in time from your sound file you will get 8 uniformly spaced frequency bins in your FFT ranging from DC up to 7/8 of your sampling rate, and each bin will have a frequency resolution of 1/T Hz where T is the time duration of your 8 samples. I think you will see from that, for your application, you will want a lot more samples (assuming you are interested in seeing more frequencies). Also learn about windowing where you taper the magnitude of your analog signal prior to taking the FFT and how that will effect the dynamic range and frequency resolution.


Cascading conventional filter with matched filter The most upvoted answer to this question suggests to use a "conventional" filter in cascade with a matched filter to first remove out of band noise (with conventional filter) and then optimize signal to noise energy (with matched filter).

Given a sampled signal where it is known that all relevant signal information is located in signal's baseband such as this one (the relevant signal here has a gaussian shape):

and where additive noise information is not known (I know it is for that specific signal, but let's assume it could have any property (colored, white, correlated, etc)).

Is it always better to use the proposed cascaded filters (conventional and then matched filter) or would it be equivalent to only use a matched filter with a gaussian shape?

No, it is not “always” better.

One would  prefilter in circumstances where it makes sense, like Eckart filtering where the signal is random and the background has a known fixed component.

Filtering, generally introduces correlation in a time series.  A simple matched filter is derived under the assumption of a white background noise. If the prefilter is a whitening filter, it is better but each situation requires its own analysis.

Food for thought, if both filters are LTI, a prefilter can be a post filter.


why Sampling with higher than nyquist frequency is causing aliasing? fm=20000 and I am sampling following signal at 65000 Hz:

x=0.5*sin(2*pi*(fm-2500)*t)+0.6*sin(2*pi*(fm-1250)*t*1.05)+0.7*sin(2*pi*fm*t*1.1)+0.8*sin(2*pi*(fm+1250)*t*1.15)+0.9*sin(2*pi*(fm+2500)*t*1.2);

and looking at frequency spectrum is not as it should look like.
Please find attachment.
It has code and result also.

The effect you see has nothing to do with aliasing. First of all, due to spectral leakage the maxima of the spectrum do not exactly reflect the amplitudes of the sinusoids. Second, note that the DFT (implemented by the FFT) only samples the actual spectrum (the discrete-time Fourier transform, DTFT) of the signal. Depending on the choice of the sampling frequency in your example, the actual spectrum (DTFT) and the location of the sample points of the DFT shift with respect to each other, and there's no guarantee that the samples of the DFT occur at (or close to) the maxima of the actual spectrum.


reduce harmonics produced by microcontroller I'm generating an 18kHz sine signal using a dsPIC microcontroller at 44100kHz (I'm not sure I can call it sampling rate, maybe generation rate, although the concept is the same).
The dsPIC's DAC then feeds an LM386 amplifier and then an 8-ohm speaker.

Since the generation method is digital, I get harmonics (36kHz) which get aliased, and I hear the parasitic low frequency (8.1kHz) which is something I wish to get rid off.

I was told that I should put an LPF with a sharp cutoff just above 18kHz at the output of the microcontroller so that those high frequency harmonics don't get amplified, but my intuition tells me that I should put an HPF, since by the time the signal is out of the microcontroller it already contains those low frequencies and an LPF will not help.

Am I correct?
Thanks!

No, you need a low pass filter aka reconstruction filter.

In the general case your filter needs to reject everything at or above Nyquist = Fs/2 = 22.05 kHz. Typically it might have have a roll-off from say -3 dB at 20kHz to -80 dB at 22 kHz. This would be quite a demanding filter to implement with analogue components however, so if you are really only interested in generating a single tone at 18 kHz then you can use a less aggressive LPF which has a -3 dB point at say 20 kHz and a stop-band which begins at 35 kHz.


Recursively obtaining all-pole coefficients from complex cepstral coefficients? This is a problem from the book Discrete-Time Speech Signal Processing by Thomas Quatieri, Chapter 6 Problem 13:

  Suppose $h[n]$ is all-pole (minimum-phase) of order $p$. Show that the all-pole (predictor) coefficients can be obtained recursively from the first $p$ complex cepstral coefficients of $h[n]$.

I was only able to solve for the opposite. i.e. complex cepstral coefficients from the pole coefficients.

If $h[n]$ is the impulse response of an all-pole minimum-phase filter of order $p$, then we have

$$\mathcal{Z}\{h[n]\}=H(z)=\frac{G}{1-\displaystyle\sum_{k=1}^pa_kz^{-k}}\tag{1}$$

where $G$ is some gain value, and $a_k$ are the predictor coefficients. Let $\hat{h}[n]$ be the complex cepstral coefficients of $h[n]$. If you know how to compute $\hat{h}[n]$ from $a_k$, then you're probably familiar with the following recursion:

$$\hat{h}[n]=a_n+\frac{1}{n}\sum_{k=1}^{n-1}ka_{n-k}\hat{h}[k],\quad 0<n\le p\tag{2}$$

From $(2)$ you can formulate a recursion for $a_k$:

$$a_n=\hat{h}[n]-\frac{1}{n}\sum_{k=1}^{n-1}ka_{n-k}\hat{h}[k],\quad 1<n\le p\tag{3}$$

with the initial value $a_1=\hat{h}[1]$.


Setting Parity Check Matrix for LDPC channel coding I am trying to implement the LDPC code for OFDM system, where I input data should be binary of length $1024$ and the output coded data of length $2048$ (means of rate 1/2).
As known in Matlab integrated function, we build the encoder and decoder, respectively, as:
ldpcEncoder = comm.LDPCEncoder;  %Encoder
ldpcDecoder = comm.LDPCDecoder;  %Decoder

After checking that, I see that parity check matrix is set by default into $32400$ x $64800$ which means the input data length must be a column of length $32400$.
I need the input data to be set manually, for example of length $1024$, how can I modify the parity check of the above encoder?
Thank you

As specified in documentation, using 'ParityCheckMatrix' you can configure the Parity Check Matrix (PCM) during the constructions of the encoder/decoder objects.
The syntax is encoder = comm.LDPCEncoder('ParityCheckMatrix',pcm) or simply
encoder = comm.LDPCEncoder(pcm);

where pcm is the desired PCM which must be sparse type.
An example for (probably poor performance) 1024x2048 PCM:
pcm = sparse([eye(1024) eye(1024)]);
encoder = comm.LDPCEncoder(pcm);
decoder = comm.LDPCDecoder(pcm);

%data
m = randi([0 1],1024,1);
%encode then decode. Note the convention of the output of comm.LDPCDecoder
y = 1 - decoder(encoder(m));
%error which equals to 0
biterr(y,m)


AM-AM plot (input power vs output power) for power amplifier I have modelled a power amplifier using weiner's model in matlab.
But when I finally plotted the input vs output power plot I got a plot that grows linear , then reaches maxima and then starts falling. As far as I know it should saturate after certain input power. Could someone tell me the reason for such behaviour? The input signal I am using is 20MHz LTE signal.

I would assume that you are only monitoring the power in-band, and what is then occurring as the signal goes into saturation is that power is being translated to other frequencies due to the non-linearity. I would suspect that this would dominantly be power in harmonic frequency locations, notably the third and fifth harmonics.


Frequency modulation in MATLAB I am trying to frequency modulate a cosine signal(for assignment). I have written the following code but not getting the required output.

dt = 0.0002;
t = 0:dt:0.05;
fm = 50;
fc = 1000;
beta1 = 1;
beta2 = 2;

%making sine and m(t)
m_t = 5.*cosd(2.*pi.*fm.*t);
sin_fm = sind(2.*pi.*fm.*t); %integral of m_t. A_m adjusted in beta

%fmmod signal
fmmod1 = cosd(2.*pi.*fc.*t + beta1.*sin_fm);
fmmod2 = cosd(2.*pi.*fc.*t + beta2.*sin_fm);

%frequency_deviation
kf_1 = (beta1*fm)/5; %5 is amplitude of message signal(Am)
kf_2 = (beta2*fm)/5;

%frequency modulated signal of matlab
fmmod1_m = fmmod(m_t,fc,1/dt,kf_1);
fmmod2_m = fmmod(m_t,fc,1/dt,kf_2);

figure(1);
plot(t,fmmod1);
title('Without using fmmod beta=3');

figure(2);
plot(t,fmmod2);
title('Without using fmmod beta=5');

figure(3);
plot(t,fmmod1_m);
title('With using fmmod beta=3');

figure(4);
plot(t,fmmod2_m);
title('With using fmmod beta=5');

I have used sind in place of sin in the code because for earlier code of Amplitude Modulation sin was not giving the required waveform but when I used sind it was giving the required waveform.

AM code

dt1 = 0.0001;
fc1 = 500*1000;
fm1 = 5*1000;
t1 = 0:dt1:0.1;

%making cosine for message and carrier for higher frequency
cos_fm1 = cosd(2*pi*fm1.*t1);
cos_fc1 = cosd(2*pi*fc1.*t1);

%making Udsb_am signal for higher frequecny
Udsb_am2h = 2.*cos_fc1 + cos_fc1.*cos_fm1;
Udsb_am1h = 1.*cos_fc1 + cos_fc1.*cos_fm1;
Udsb_am_half_h = (0.5).*cos_fc1 + cos_fc1.*cos_fm1;

figure(1);
plot(t1,Udsb_am2h);
title('AM for Ac=2');
figure(2);
plot(t1,Udsb_am1h);
title('AM for Ac=1');
figure(3);
plot(t1,Udsb_am_half_h);
title('AM for Ac=0.5');

What I want to ask from the question?
When should I use sind and sin?(I know sind is for degree and sin is for radian).Asking question about sind and sin because when I saw the code on internet they have used sin not sind.

What is wrong in the code of FM?

The difference is the expected input. The cosd function expects the input to be expressed in degrees, and the cos function expects the input to be in radians. So you have cosd(theta)==cos(deg2rad(theta)). This is just from the MATLAB documentation page.

Are you using these correctly? For example, you do cosd(2*pi*F*t). So we have $\frac{2\pi \text{ radians}}{\text{cycle}}.\frac{F \text{ cycles}}{\text{second}}.[\text{seconds}]=2\pi F \text{ radians}$, as the input. You should be using cos not cosd because your input is expressed in radians.

BTW, I do not check the correctness of the actual code. That is on you.


Upsampling vs downsampling. Which to use when? Downsampling reduces dimensionality of the features while losing some information. It saves computation.
Upsampling brings back the resolution to the resolution of previous layer.
My question is which is better or when to use which one?

You should use the one you need for your problem, when you know which components of your signal are of interest to you.
Let's say you have in your electronic editing an ADC digitizing 40M samples per second to study a heart rate of 70 beats per minute, you are very likely to work with useless information, that's why it will be better to down-sample your signal.
In another case, it can be very useful to oversample, some algorithms to compute an FFT work well with sample sets of power 2, it may be necessary to oversample your signal to apply a fast FFT.
That is why you should know the characteristics of your signal of interest before applying any of these functions.


How to Normalize Frequency in MATLAB I know all about conversion of frequency from HZ to normalized frequency and from normalized frequency to angular frequency ,but I wonder how it is done in MATLAB.
I have the following code:

f:0:0.001:1;
omega=pi*f

here f=normalized frequency and omega=angular frequency. My question is how do we get the interval value(0.001) for our normalized frequncy.

Thanks for any help

The interval resolution should be determined according to data you want to analyze.

It is equivalent to the question, what magnification of your Micro Scope would you use?
Well, it depends on the size of the things I want to look at.

If your data within the sampling frequency you set has different energy in each frequency bin and you'd like to see it, use small intervals.

Though,
The real question, how do you intend to use the Frequency Vector?


Delay and sum Beamforming What is the difference between the weight vector (w) and the steering vector a(θ)?

I initially assumed they were the same thing, but the literature seems to be to be a bit fuzzy with this (or maybe its just me).

They're essentially the same, with the steering vector typically representing phase shifts (==delays) using complex numbers from the unity circle, and the weight vector often containing real values to scale the samples. In practice, both can be unified into one and hence, the terms are sometimes used interchangeably. Also note that for many algorithms, the individual amplitude of the samples doesn't matter, and hence, the steering vector's magnitudes can be chosen arbitrarily (as long as non-zero).


Finding out modulation index and DC offset I have a question form my teachers, and I cannot understand why I can find out the modulation index form the figure.

The question provide a Figure like this:

  And the information signal is a sinusoidal test signal with peak amplitude $6\ \rm V$ and is applied to an AM-DSB-C modulator, the Fourier spectrum of the modulated signal is shown above.

The solution is like this:

  As $\displaystyle \frac {\frac {Acm}{4}}{\frac{Ac}{2}}=\frac{3}{4}$, $m=1.5$, so the modulation index is $1.5$
  Moreover as $\displaystyle  m=\frac{x}{c}$, the peak amplitude of $s(t)$ is $6\ \rm V$, so the DC offset($c$) is $4\ \rm V$.

I know where is the $\displaystyle \frac {Ac}{2}=4$ come form

as $s_{AM-DSB-C}(t)=A\big(s(t)+c\big)\cos(2 \pi f_c t)$

\begin{align}
\mathcal F\left\{s_{AM-DSB-C}(t)\right\} &=S_{AM-DSB-C}(f)\\
&=\frac{A}{2}\big[S(f-f_c)+S(f+f_c)\big] + \frac {Ac}{2} \big[\delta (f-f_c) + \delta (f+f_c)\big]
\end{align}

so, form the second term we can get $\displaystyle\frac {Ac}{2}=4$

But, I cannot understand how the solution can get $\displaystyle \frac {Acm}{4}=3$.

The sinusoid test signal is given by
$$
s(t) = V_0\cos(2\pi f_0)
$$
and its Fourier transform by
$$
S(f) = \frac{V_0}{2}\big[\delta(f-f_0) + \delta(f+f_0)\big].
$$
Inserted in $S_\mathrm{AM-DSB-C}$ it yields
$$
S_\mathrm{AM-DSB-C}(f) = \frac{AV_0}{4}\big[\delta(f-f_c-f_0) + \delta(f-f_c+f_0) + \delta(f+f_c-f_0) + \delta(f+f_c+f_0) \big] + \frac {Ac}{2} \big[\delta (f-f_c) + \delta (f+f_c)\big].
$$
Comparing that with the figure given we see that $AV_0/4 = 3$ and $Ac/2 = 4$. Additionaly, $f_0 = 10\text{ kHz, } f_c = 100\text{ kHz}$. Using $m = V_0/c$ and $V_0 = 6\text{ V}$ we find
$$
m = 1.5, c = \frac{V_0}{m}=4\text{ V}.
$$


What is the difference between an equalizer and a channel estimator? I was wondering what is the difference between channel estimator and channel equalizer?
Thanks!

In communication systems, transmitted signals are distorted by the physical medium (the channel) charactheristics.

The channel estimator tries to identify the transmission channel characteristics, by computing (estimating) various of its model parameters.

Then, the channel equalizer tries to restore to distorted (received) signals back into their originals via some sort of inverse filtering (equalization) on the received signals, which was designed according to the estimated channel parameters.

Both operations are approximate; depending on the application, they can be very effective or not so.


Understanding LPC for Formant Estimation I went through the Matlab tutorial on Formant Estimation using LPC Coefficients. Though I vaguely understand the details, it's not entirely clear why we need to do this. From http://person2.sol.lu.se/SidneyWood/praate/whatform.html:

  A formant is a concentration of acoustic energy around a particular frequency in the speech wave

Why is it not enough to take the DFT of the audio signal (after some pre-processing if necessary)? In the frequency-domain, the peaks correspond to these concentrations, correct?

The frequency resolution of DFT is limited to the number of time samples. On the other hand proper LPC can have high resolution.