texts: How is `aicbic` or other matlab function used to find the AR order of specific data?

texts:
let us consider following matlab code

help aicbic
 aicbic Akaike and Bayesian information criteria

  Syntax:

    [aic,bic] = aicbic(logL,numParam,numObs)

  Description:

    Given optimized log-likelihood function values logL obtained by fitting
    a model to data, compute the Akaike (AIC) and Bayesian (BIC)
    information criteria. Since information criteria penalize models with
    additional parameters, aic and bic select models based on both goodness
    of fit and parsimony. When using either AIC or BIC, models that
    minimize the criteria are preferred.

  Input Arguments:

    logL - Vector of optimized log-likelihood objective function values
      associated with parameter estimates of various models.

    numParam - Number of estimated parameters associated with each value
      in logL. numParam may be a scalar applied to all values in logL, or a
      vector the same length as logL. All elements of numParam must be
      positive integers.

  Optional Input Argument:

    numObs - Sample sizes of the observed data associated with each
      value of logL. numObs is required for computing BIC, but not AIC.
      numObs may be a scalar applied to all values in logL, or a vector the
      same length as logL. All elements numObs must be positive integers.

  Output Arguments:

    aic - Vector of AIC statistics associated with each logL objective
      function value. The AIC statistic is defined as:

        aic = -2*logL + 2*numParam

    bic - Vector of BIC statistics associated with each logL objective
      function value. The BIC statistic is defined as:

        bic = -2*logL + numParam*log(numObs)

Let us suppose that we have the data in a matlab vector X.

How do I  apply this matlab function (aicbic) to X?

I am trying to find the AR order of the data.


texts:
You are testing a new device to determine its characteristics. You find the magnitude of its frequency response shows a peak at 500 Hz, and the measured overshoot of the step response is 40% of the steady state level. Find the system’s characteristics and confirm your answers by generate plots of its frequency and step response using either Simulink or Pspice.

I'm trying to figure out how to find the natural frequency with the given info.

$\frac{500\,\text{Hz}}{140\,\%}\approx 357.145$;  I'm assuming my $r_o = 357.145$, where we settle at.

Given percentage overshoot (40%) I can solve for our damping ratio $\zeta = 0.279$.   I have formulas to solve for Natural Frequency ($\omega_n$) but not enough known variables to solve for it.  Does anyone know how to solve for $\omega_n$ here?

In our case, $\beta = 1.288$, if that helps.


texts: Is an FIR filters 'transition width' frequency dependent?

texts:
I was recently brushing up on some FIR filter design and was browsing the Wikipedia page. A claim in the article caught my attention:

  The main disadvantage of FIR filters is that considerably more
  computation power in a general purpose processor is required compared
  to an IIR filter with similar sharpness or selectivity, especially
  when low frequency (relative to the sample rate) cutoffs are needed.
  However many digital signal processors provide specialized hardware
  features to make FIR filters approximately as efficient as IIR for
  many applications.

Article

Now what caught my attention was the statement made about how more computation is required when a low cutoff frequency relative to the sample rate (i.e. low normalized cutoff frequency) is specified.

My question: is this true? Or is it just a poorly written statement that assumes you need sharp transitions when working with low cutoff frequencies?

I was always been under the impression that FIR filters become very difficult to implement when trying to achieve sharp transitions, not certain frequencies. Have I been wrong in this assumption?

Furthermore, this article seems to confirm the fact that the decay rate of the impulse response is independent of a chosen cutoff frequency. So assuming the impulse response has a constant decay rate, shouldn't an FIR filter of fixed length being equally good (or bad) regardless of cutoff frequency?


texts: Obtain SIB1 messages from CPRI the captured DL IQ data

texts:
I am a beginner in the signal processing, and I know some python and I have a project I'm working on which is about Obtaining SIB1 messages from CPRI the captured DL IQ data.

I have captured some DL IQ data from CPRI Prob between the LTE BBU & RRU, and I'm interested in decoding the messages and get some basic info, like band information, or band indicator, just to know which band this AxC is related to.

For example, If the CPRI link has 8 AxC, antenna carriers, for two bands, B1 has 4 AxC, and the other second band has the remaining 4 AxCs.==> I'm interested in decoding the SIB1 messages from this captured DL IQ data so I can know which AxC is belong to which band.

The captured data coming in two csv files: (each band in one seperate file).
Both files contain IQ data from 4 antennas on each band. Each column in a file belongs to one antenna.
There are 308736 samples for each antenna
Sample rate is 15.36 Megasamples/second.

From my research today, I found some provided answers recommending GNU Radio for such project, so it will be appreciated if any one of you can provide a feedback, advise, or recommendation, and make my day easier:

How to go from IQ to SIB1 messages.
Is GNU Radio simulation is the best to be learned/used or some other python packages can help.
At the end I want my CPRI prob know this automatically, instead of manual download the csv capture and analyse to see which AxC antenna related to which band.

One document I found about GNU Radio simulation: A LTE Receiver Framework Implementation in
GNU Radio


texts: How to compute autocorrelation of signal defined by difference equations?

texts:
I have no experience with difference equations and I want to learn how to compute the following, but I found no resource online. Any help would be greatly appreciated.

Find:

$$\mathbb{E}\left[d[n]d[n + k]\right]$$

Where:

$d[n]$ is a discrete time signal of the form:

$$d[n] = \alpha \,d[n−1]+v[n] \quad (\alpha \in [0;1])$$

$v[n]$ is a sequence of uncorrelated zero-mean unit-variance random variables.


texts: Spacing between gaussian windows for STFT

texts:
I'm computing discrete short time Fourier transform. Data is split into overlapping chunks and gaussian window is used for each chunk. However, I'm not sure how much overlap there should be between chunks.

If they are too sparce I will loose information between them (gaussian window for each chunks is shown here):

On the other side, pack them to densely and I will do a lot of redundant computations:

Something in between is necessary:

But what exactly should it be?


texts: Can FFT be used to find peaks and pits on audio files

texts:
I'm able to read wav files and its values. I need to find peaks and pits positions and their values.

First time, I tried to smooth it by $\frac{(i-1 + i + i + 1)}{3} $
formula, then searching on array as array[i-1] > array[i] & direction == 'up' --> pits style solution. But because of noise and other reasons of future calculations of project, I'm trying to find a better working area.

Since couple days, I'm researching FFT. As my understanding, FFT translates the audio files to series of sines and cosines. After FFT operation the given values is a0's and a1's for $a_0 + a_k * cos(k*x) + b_k * sin(k*x)$ which k++ and x++ as this picture:

My question is, does FFT helps to me find peaks and pits on audio? Does anybody has a experience for this kind of problems?


texts: What Approximation Techniques Exist for Computing the Square Root?

texts:
I have very limited resources as I'm working with a microcontroller.  Is there a taylor-series expansion, common lookup table, or recursive approach?

I'd prefer to do something without using math.h's sqrt()

http://www.cplusplus.com/reference/cmath/sqrt/


texts: What is the magnitude of the impulse function at t=0?

texts:
The unit impulse function, in a few textbooks that I have referred, has a value of 0 at t≠0 , and an area of unity (1). The height of the impulse function also tends to infinity at t=0. But since it's height cannot be it's magnitude (as we generally use it's unit area in calculations), what would be the value/amplitude/magnitude of the unit impulse function at t=0?

Before you say that my question has answered itself, and that it's value is one(1) at t=0, note that this is not defined as such in any standard textbooks. This question has been bothering me for quite a while now, so any attempt at elaborating a solution will be highly appreciated.


texts: How can I smoothly interpolate between 2 position?

texts:
I've got a 1D signal (position of a servo motor over time) and I've extracted 'peaks'/'key' positions picking running average "local extrema" points.

Below is are 2 plots from 2 servos and the white markers indicate key positions I'd like to interpolate between:

I've played with linear interpolation, but the motion looks robotic. I'd like to get a more natural/smoother motion and I imagine it's posible by using bezier interpolation between two key points. These two points would be the start and end point of the bezier I imagine, but I need to solve/find out the 'anchor' points as well. Since I have all the points I can get the velocities as well I imagine, but not sure how to plug them into the equation.

Any idea on how I can find the other two members of a bezier equation for easing/smooth interpolation ?

Also, since I'm new to this: is my approach ok ?  Is there a simpler solution I'm missing ?


texts:
I am trying to determine the synchronicity between two oscillating signals. I was attempting to do this by making use of the cross correlation between the two signals and looking at the maximum value. Is this a correct approach?

If so, which would be the better approach, to normalize the signals between 0 and 1 or to normalize the cross correlation?

I don't really mind about the difference in amplitude but I would like to see levels of synchronicity. Originally I thought to do cross correlation and look at the max of the graph however, I am not sure that this is the correct approach. I have also been looking at coherence of the spectra


texts:
I was looking at a technology which enables virtual advertisement on the perimeter ads in soccer. It's relatively an old technology for a few years back, but I'm surprised to see the quality of it. They overlay the physical perimeter ads with new ads on different camera feeds. Basically from this link. And here is a video on Youtube.
LED perimeter system with integrated virtual replacement technology - makes it possible to change the content seen by different viewing audiences around the world
How does this overlaying mechanism work? Is there any hardware involved (like the LED boards itself), or it is purely based on Computer Vision and Augmented Reality? How can it be such high-quality without glitches purely based on computer vision?!


texts:
I am attempting to design a low-pass filter for a subwoofer, as the audio amplifier that is connected to it does not filter any of the sound coming in, and it is getting pumped full of high frequency sound. I have settled on a Butterworth LC filter, and want a crossover frequency of around 200Hz. I have determined that I want a capacitor somewhere in the range of 140-160 uF and an inductor in the range of 4.5-6.5 mF to hit that frequency, and from a circuit standpoint, it seems to make sense. Source: this website

Here is where I hit a difficulty. I am aware that I don't want an inductor with very much DC resistance and is rated for a high current throughput (possibly in excess of 5A). I don't know which to choose.

I also am at a loss for the capacitor to choose, as I need one which is two-way (able to handle AC), but I am not sure for the other requirements of it.
If anyone can clarify these points for me, and give suggestions as to the actual parts to use, I would be forever grateful. It is also worth noting that I am on somewhat of a budget, so please try to recommend affordable products, even if they have lower performance.


texts: Seeking compressive sensing imaging demo in MATLAB

texts:
I am looking to experiment in MATLAB with image restoration/reconstruction using compressed sensing. I am relatively new to compressive sensing, and was inspired by the following Wired article to learn more about it:

Fill in the Blanks: Using Math to Turn Lo-Res Datasets Into Hi-Res Samples

In the Wired article linked above, the authors show a low-res picture of US President Obama, and snapshots of the image as it progresses through various stages of "recovery," ultimately ending up with a clean-looking image.

I have read quite a few papers, and even have coded a simple algorithm for sparse signal recovery in MATLAB, but I would like to be able to "play around" with a low-quality image in MATLAB, using some compressed sensing algorithms. I think it would help solidify my understanding of what is going on by incorporating a visual example.

So does anyone know of any such MATLAB demos/examples which provide the MATLAB code and a sample image (or images) that would permit one to experiment with the recovery process?


texts: Time-differentiation property of Fourier transform for $\sin(\omega_0t)$ and $\cos(\omega_0t)$

texts:
As I know, the Fourier transform has the below property which is called time-differentiation:

$$
\frac{dx(t)}{dt}\leftrightarrow j\omega X(j\omega)
$$

and the Fourier transform of the cosine and the sine signals are $\pi(\delta(\omega-\omega_0)+ \delta(\omega+\omega_0))$ and $\frac {\pi} {j}(\delta(\omega-\omega_0)- \delta(\omega+\omega_0))$, respectively.

The derivative of $\sin(t)$ is $\cos(t)$. However, the Fourier transform of $\cos(t)$ cannot be obtained from multiplying the Fourier transform of $\sin(t)$ with $j\omega$.

could any one help me.


texts:
A signal $$x(t)=\left|\sin(2\pi(1000)t)\right|$$ ($\left|x\right|$ is the absolute value function) is fed to an ideal low-pass filter with cutoff frequency of $1500\,\mathrm{Hz}$ to produce an output $y(t)$, what is the RMS voltage of $y(t)?$

My initial belief is that it will be zero since the absolute value doubles the frequency of the signal and so the frequency of $y(t)$ will be $2000\,\mathrm{Hz}$ which is beyond the cut-off frequency of the ideal low pass filter.

However, upon reflecting on the Fourier transform of $y(t)$, I believe it will have a DC component with zero frequency which will give it an RMS voltage. Moreover, $y(t)$ is a fully rectified signal which, when passed through low pass filters in reality (i.e. capacitors) will yield a smooth DC signal.

How do I find the RMS voltage of the resulting $y(t)?$ Is it really nonzero? Will it be just the same as the RMS voltage of $x(t)?$


texts:
I have a simple question, that might not be up to the level in this forum:

  In theory, a signal $s(t)$ that is received through a multipath channel reads:
  $$r(t) = \sum\limits_{k=1}^p \alpha_ks(t - \tau_k)$$
  where $r(t)$ is the received signal,  $\alpha_k$ is the gain of the $k^{th}$ path and $\tau_k$ is the delay of the $k^{th}$ path.

But

  In MATLAB, we have $s[n]$, sampled at some sampling period, $T_s$. To simulate $r[n]$, the only problem faced is that $\tau_k$ will not be a multiple of $T_s$. My question is "Is it a good idea to upsample the signal $s[n]$ up to a desired sampling period, say $NT_s$". If so, why should filtering be done after the upsampling process ?


texts: Understanding the FFT phase spectrum with a simple example

texts:
I'm trying to compute the DFT using scipy's functions. I don't understand why the phase spectrum of a simple sine wave with 2 Hz frequency doesn't show $\pm\pi/2$ at the $\pm 2Hz$ frequencies. Instead, the phase plot seems to have some linear dependence in the frequency, which I don't understand. I provide the code for assistance. How can this be fixed? Looks like a simple issue I'm not grasping. Please help.

import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import fft, fftfreq

# 1Hz sine wave
npts  = 100
tmax  = 10
t     = np.linspace(0, tmax, npts)
y     = np.sin(2*np.pi*2*t)
dt    = tmax/npts

# FFT computation
Y      = fft(y)
freq   = fftfreq(npts, d=dt)
amplit = abs(Y)/npts
phase  = np.angle(Y)
phase  = phase / np.pi

fig, ax = plt.subplots(1, 2, figsize=(10,4))
ax[0].plot(freq, amplit)
ax[1].plot(freq, phase)

ax[0].plot([2,2], [0,0.5], '--')
ax[0].plot([-2,-2], [0,0.5], '--')
ax[1].plot([2,2], [-1,1], '--')
ax[1].plot([-2,-2], [-1,1], '--')


texts: Suitable sampling rate for triangle wave

texts:
I know this is a fairly simple question but I can't convince myself of the right answer.

For a 5KHz triangle wave, being sample at Fs sampl/second, what's a suitable choice for Fs such that an accurate reconstruction of the cont. time wave can be generated?

My first thought was to satisfy Nyquist, so sample at twice at the wave rate. But, is there something I'm missing here? Like for the square wave, where technically the sample rate is infinite for a perfect reconstruction, is there something similar here?

Thanks!


texts:
I just need some advice on a project. I have an unwanted sinusoid that has time-varying frequency (its not sinusoidally varying) which is relatively constant over a small section. I've been advised to use a Linear Predictor in order to deal with it and remove it.

I just need some advice on how to go about implementing it. I understand the theory, my main issue is with this equation for a forward-linear predictor

$$
\hat{x}(n) = \sum_{i=1}^p a_i\,x(n-i)
$$

Now, what's the point in estimating the next value of the signal if it's also going to contain the same unwanted sinusoid component? I understand how to use this (finding coefficients using Yule-Walker/Levinson-Durbin)

EDIT: More info, Sinusoid is added on top. I can almost hear the telephone conversation clearly when the sinusoid dips down. I understand that for a small enough section, I could just use a notch filter but of course it's not feasible to do this for each small section manually. I need a way of predicitng it's frequency hence the linear predictor The sinusoid has a significantly higher amplitude than the conversation and it's frequency varies from 200Hz to 2kHz


texts:
I have been reading a couple of definitions on signal such as:

  Any quantity measurable through time over space or any higher
  dimension can be taken as a signal.

And this,

  In the physical world, any quantity exhibiting variation in time or
  variation in space (such as an image) is potentially a signal that
  might provide information on the status of a physical system

Taking human voice as example,

A sound signal is converted to an electrical signal by a microphone, generating a voltage signal as an analog of the sound signal, making the sound signal available for further signal processing. Sound signals can be sampled at a discrete set of time points.

I can understand where times comes in but where's the "space" come here according to the definition of signal?


texts: Lightweight test for bimodal distributions?

texts:
Here, and in the stats stackexchange, seem to be answers that reference tests for bimodal distributions that involve iterative binning or iterative curve fitting methods.  However "eyeballing" a plot of a data set often shows a clear bimodality (say a 10 dB dip or several standard deviations  between two clear mode peaks, etc.), versus a single "hump", or something ambiguous (less than a 3 dB dip).

Are there any lightweight algorithms (computationally efficient, single pass, or deterministic low iteration count tests) which can separate out the easy cases (clearly unimodal, or clearly bimodal), versus data sets that may require more computationally intensive tests for likelihood of multi or bi-modality?

Or, in signal processing terms, how can one non-graphically and quickly test whether some area in a noisy spectrogram is more likely to contain 2-tone FSK (with an unknown separation), rather than noise or some other kind of modulated signal?


texts: Bilinear interpolation implemented by convolution

texts:
I read the paper Deep Feature Flow for Video Recognition https://arxiv.org/abs/1611.07715.

In Sec.3, the author implements bilinear interpolation like this:

$$f_i^c(p)=\sum\limits_{q}G(q,p+\delta p)f_k^c(q) \tag{1}$$

Where $q$ is the point from the source image, and $p$ is the points on the target image. $\delta p$ is the distance the point moved each point $p$ (not $\delta \bullet p$). $G$ is defined as

$$G(q,p+\delta p)=g(q_x,p_x+\delta p_x)g(q_y,p_y+\delta p_y)\tag{2}$$

And the bilinear interpolation is defined in wiki as:

$$f(x,y)\approx {\frac {y_{2}-y}{y_{2}-y_{1}}}f(x,y_{1})+{\frac {y-y_{1}}{y_{2}-y_{1}}}f(x,y_{2})\tag{3}$$

I think the operation $(1)$ and $(3)$ is equivalent. How can I derive the filter $(1)$ from $(3)$?


texts: getting different output when using FFTW3 Library than using Rosetta FFT

texts:
I have the following array containing sine waveform that I want to compute FFT:
{ 0.0, 0.114124, 0.370557, 0.0, -0.576201, -0.370557, 0.0, 0.0 };

When using FFTW3 Library and running those 2 functions:
plan_forward = fftw_plan_dft_r2c_1d ( N, in, out, FFTW_ESTIMATE );
fftw_execute ( plan_forward );

I get the following output:
-0.462077
0.918921
-0.946758
0.233480
0.050790
0.000000
0.000000
0.000000

When using Rosetta implementation to FFT I get the following output:
-0.462077
0.918921
-0.946758
0.233480
0.050790
0.233480
-0.946758
0.918921

As you can see, some of the output is simillar between the 2s and some are different, I can't figure out why.
Here is my full code:
#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <fftw3.h>
#include <iostream>
#include <cmath>
#include <fstream>
#include <complex>
#include <iostream>
#include <valarray>
using namespace std;

const double PI = 3.141592653589793238460;

typedef std::complex<double> Complex;
typedef std::valarray<Complex> CArray;

// Cooley–Tukey FFT (in-place, divide-and-conquer)
// Higher memory requirements and redundancy although more intuitive
void fft(CArray& x)
{
    const size_t N = x.size();
    if (N <= 1) return;

    // divide
    CArray even = x[std::slice(0, N/2, 2)];
    CArray  odd = x[std::slice(1, N/2, 2)];

    // conquer
    fft(even);
    fft(odd);

    // combine
    for (size_t k = 0; k < N/2; ++k)
    {
        Complex t = std::polar(1.0, -2 * PI * k / N) * odd[k];
        x[k    ] = even[k] + t;
        x[k+N/2] = even[k] - t;
    }
}

int main(void)
{
    int i;
    double y;
    int N=8;
    double Fs=30000;//sampling frequency
    double  T=1/Fs;//sample time
    double f=5000;//frequency
    double *in;
    fftw_complex *out;
    double t[N-1];//time vector
    fftw_plan plan_forward;

    in = (double*) fftw_malloc(sizeof(double) * N);
    out = (fftw_complex*) fftw_malloc(sizeof(fftw_complex) * N);
    for (int i=0; i< N;i++)
    {
        t[i]=i*T;
        in[i] =0.7 *sin(2*M_PI*f*t[i]);// generate sine waveform
        double multiplier = 0.5 * (1 - cos(2*M_PI*i/(N-1)));//Hanning Window
        in[i] = multiplier * in[i];
    }

    Complex test[] = { 0.0, 0.114124, 0.370557, 0.0, -0.576201, -0.370557, 0.0, 0.0 };
    //const Complex test[] = (Complex*)malloc(sizeof(Complex) * N);
    //CArray data(test, 8);

    printf ( "\n" );
    printf ( "  Input Data:\n" );
    printf ( "\n" );

    for ( i = 0; i < N; i++ )
    {
        printf ( "  %4d  %12f\n", i, in[i] );
        test[i].real(in[i]);
    }
    CArray data(test, 8);

    fft(data);
    plan_forward = fftw_plan_dft_r2c_1d ( N, in, out, FFTW_ESTIMATE );
    fftw_execute ( plan_forward );

    printf ( "  output Data from FFTW3 Library:             output Data from my FFT implementation\n" );
    for ( i = 0; i < N; i++ )
    {
        printf ( "  %4d  %12f           %12f\n", i, *out[i], data[i].real());
    }
    return 0;
}


texts:
I'm using a Raspberry Pi 2 and an AD/DA Board, with no extra circuit. Reading on the internet though, I found that DACs need a low pass filter in their output, and I don't understand if the behaviour of my DAC is normal, and I wonder if I need the low pass filter.
Its output goes down to 0V at every sample, is it normal? Is this the reason why low pass filters are used?

Thanks in advance,

Alberto

I tried with a better oscilloscope to analyze the signal, and it's all fine. I'm so sorry if I created a useless thread, but I still have learnt new terms I didn't know before.


texts: Objective, repeatable sound quality measurements and metrics

texts:
I am currently doing an electronics design for a basic 2.5W speaker amplifier for an embedded system based on COTS parts (e.g. codec and amp), similar to what one would find in a tablet.

Aside from a wide range of parts to choose from, there are other engineering trade-offs, such as the filter on the CLass D amp for EMI/EMC and the noise on the power supply rails.

This, of course, is even before deciding on a speaker.

Is there a standard method for an objective and repeatable evaluation process? If so, what sort of equipment and methods do I need, and what metrics do I need to calculate, and how to these metrics relate to the perception of audio quality?

(Yes, of course, there is SNR THD and frequency response, but I am yet to relate these to how a human perceives the quality of the audio).

(Second note: This isn't an evaluation of codecs, compression, sampling or other DSP artefacts.  I have control of the bitstream I send to the audio codec chip. Advice of choosing the signal I should send is most welcome).


texts: How to interpret these different Fourier analysis of this audio signal?

texts:
This is my first dive in DSP. I would like to familiarize myself with frequency analysis. I have two audio tracks which should be digitized at 16bit-44.1kHz and 24bit-192kHz (music, presented as a 24bit-192kHz sample) respectively.

I wanted to identify the effect of the low-pass filter around the Nyquist frequency (22.05kHz and 96kHz respectively).

Edit: I completely reworked the question.

Software used:

I basically estimated the power spectral density using Welch's method
as implemented by scipy.signal.welch in the Scipy library of the Python programming language.

Basically, I used a script equivalent to:

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

# Load data from one channel (#0) for each sample file
wav192 = readwav("24b-192khz.wav")[:,0]
wav44 = readwav("16b-44khz.wav")[:,0]

# DoE: 2 sample size and two windows types
chunks = [256, 4096]
windows = ["hanning", "boxcar"] # boxcar is rectangular

# Prepare a figure
plt.figure()
# Calculate density spectra and plot
for N in chunks:
    for w in windows:
        f, Pxx44 = signal.welch(wav44, fs=44100, window=w, nperseg=N, nfft=2*N, scaling="density")
        plt.semilogy(f, Pxx44)
plt.legend(["chunk=%d; window=%s"%(c, w) for c in chunks for w in windows])
plt.xlabel("Frequency (Hz)")
plt.ylabel("Density (I$^2$/Hz)")

The power spectral density of the 44.1kHz audio sample:

Which is basically just as expected:

The chunk size, i.e. the number of samples per fft-transform segment in the real domain, does not change the density a lot if a Hann window is used.
The chunk size effect is clearly visible with the boxcar (rectangular) window. From what I understand, this is because of spectral leakage which diminishes as the chunk size increases. Is that correct?
The low-pass filter effect at the Nyquist frequency (22.05kHz)

So far, so good.

The power spectral density of the 192kHz audio sample:

Good point:

Same behaviour in regard to the chunk size and the window. Is spectral leakage really that strong? That's pretty impressive.

Oddities:

What the heck is happening?
Where is the low-pass filter near the Nyquist frequency?
Why are very-high frequencies even increasing? Could that be related to the choice of the windowing function?

From my interpretation, there is no low-pass filter visible because basically no audio system would go above 192kHz and generally, the software/hardware creator are smart enough to apply a low-pass filter designed with regard to the actual output bandwidth of the audio system.

As for the increasing audio signal above 57kHz, I really can't explain it: the original audio sample is some classical music. I wouldn't expect any instrument to generate louder sounds in that range or frequencies. Any idea? Could this be an example of upsampling?


texts: What could be causing these humps every 10 Hz on my PSD estimate data?

texts:
I am working with EEG recordings about 1,000,000 data pts long, recorded at 4,000 Hz.

I'm generating PSD estimates for these recordings using the matlab "periodogram" function like so:

[power,f]=periodogram(testdata,hamming(length(testdata)),[],4000);

and plotted like so:

scatter(f,log(power),1)

When I look at my plotted data (frequency vs. power), I see these strange, camel-like humps with maxima at every 10 Hz. The humps kind of resemble the plot of the absolute value of a sine wave. These humps accompany and are quite distinct from what I believe to the the non-artifactual baseline data.

When I run the periodogram function on smaller portions of the data as opposed to the full 1,000,000 point recording, the humps and the baseline appear to merge together.

Does anyone know what these humps might represent?


texts: Taylor series expansion in mean shift tracking

texts:
I am trying to understand the first order Taylor series expansion when deriving the math behind mean shift tracking. Given the Battacharya coefficient between the target descriptor and candidate descriptor $\rho[\textbf{p}(y),\textbf{q}]$, how does
$$\rho[\textbf{p}(y),\textbf{q}] \approx \frac{1}{2}\sum_m\sqrt{p_m(y_0)q_m} + \frac{1}{2}\sum_mp_m(y)\sqrt{\frac{q_m}{p_m(y_0)}}$$
I am trying to see how each term can be related to the taylor series expansion about a point $a$ given by $$f(x) \approx f(a) + f'(a)(x-a)$$
What does $f$, $x$ and $a$ correspond to in the mean shift problem ?
Would appreciate any tips and pointers regarding this. Many of the lectures online do not really explain the details.


texts:
It can be shown that the image derivative in the $x$ direction is given by:
$$
\frac{\partial f}{\partial x}=\frac{2\pi i}{N} \mathcal F^{-1}\left(u\cdot \mathcal F(f(x,y)\right)
$$
where $N$ is the width of the image in the $x$ direction and $u$ denotes a variable in the fourier space ($F(u,v)$).

Now, in the image space, derivation is a convolution with this matrix:
$$\partial_x =\left(\begin{array}{cc}
-1 & 1\end{array}\right)$$

this means (from to convolution theorem) that the fourier transform of $\partial_x$ is a matrix $D(u,v)$ such that $D(u,v)=u$.

so I tried to see this in python, but something seems to be wrong here:

d = np.zeros((21,21))
d[10,10] = 1
d[10,9] = -1
D = np.fft.fft2(d)
plt.matshow(np.fft.fftshift(D))

Can anyonw see where I am wrong?


texts: Coherence bandwidth/delay spread in OFDM

texts:
Say we want to design an OFDM system and we have some parameters given to us, i.e. the coherence time $T_c$ and the max delay spread $T_m$.

It's clear to me why we need to ensure $T_s(N+N_{cp}) \ll T_c$ where $N_{cp}$ is the cyclic prefix in samples - we want the channel to stay coherent during the entire symbol. What is not clear is what happens to the constraint $T_s(N+N_{cp}) \gg T_m \Leftrightarrow B_N \ll B_c$. Apparently, in OFDM it suffices with $T_s N_{cp} \geq T_m$.

It makes intuitive sense in a way, but I'm still not sure exactly how the coherence bandwidth interacts with the subchannel bandwidths. Why don't we need to ensure $B_N \ll B_c$?


texts: How to transform a signal to go through specific points?

texts:
I have a 1d signal obtained using a Fourier based resample method (TDIFDZP) for which the resampled points don't necessarily go through the original samples. I want to transform the upsampled signal to go through these points to reduce the error between the resampled and the original signal as shown in the image below. I'm using scipy.sigal.resample implementation in python.
I'm looking for a transformation of the resampled signal so that it goes through the original samples.

Edit: A little of background, I'm exploring different interpolation methods for sub-Nyquist sampling. My original signal's Nyquist rate is ~20Hz but the number of samples available in my application is <6.


texts:
I'm trying to understand how a dynamic convolution model relates to something like a Volterra series, and what kinds of effects the latter can capture that the former can't (and vice versa). By dynamic convolution, I mean a system where the impulse response varies with the amplitude of an impulse.

I care only about dynamic convolutions that are well-behaved and "smooth" in some sense, i.e. where the impulse response varies continuously or smoothly as a function of the amplitude.

I have this intuition that a Volterra series with nonzero coefficients only on the diagonals are somehow related to dynamic convolutions. In other words, a series that reduces to the following:

$y = h_0 + h_1 \ast x + h_2 \ast x^2 + h_3 \ast x^3 + ...$

where there are never any terms like $x[t] \cdot x[t-1]$ and so on.

I'm curious if I'm on the right track with that, and in general, whether any dynamic convolution can be transformed into a Volterra series in that form or vice versa.


texts:
I'm new to the whole DSP area, but I've built a program which plays music and does FFT to produce some cool visuals for the music in realtime.

The FFT algorithm I chose was from a C# library called NAudio. I applied a hamming window before shipping it off to the FFT, from which I get an array of complex. Nothing special. Works fine.

But now I want to display this data to the user, so I converted the bins to magintudes following this formula: magnitudes[i] = sqr(real[i]² + imaginary[i]²)

This with some amplification yields a pretty good result to display for the user, but I feel that most of the levels being displayed are just random noise, I want to show the user what he/she "hears", where the melody or vocal are being amplified OVER the other background music.

The thing I'm unsure of is how to approach this problem. What can I do to filter out everything else and only have the melody/vocals left. I can understand there's not 100% clean way you can do this, but what method could give me a decent result?

Thank you


texts: Thresholding a difference image with a very high kurtosis

texts: I have an difference image produced by substracting two images of a certain location taken a two different time steps. The resulting image of this operation shows a histogram with a very large peak centered a 0. Since this naturally can be identified with no-change between dates I interpret this as a difference image showing a very small amount of change. Yet by the magnitudes of the observations located at the tails; strong change. I would like to threshold this image to identify three classes: negative change, no-change and positive change. I have been experimenting with classical thresholding techniques (Otsus, Kapurs, Rosins, etc) applied to each of the sides of this unimodal distribution to achieve this. I beleive I am not getting to good results because these techniques are not designed for a very spikey "distribution" with very long tails. Can anyone propose a method that can be suitable for this problem or some sort of preprocessing of the original data (difference image) which can help me out on this?


texts:
I have two FFTs that represent two impulse responses. I want to interpolate between the two of them to generate a 3rd FFT in the manner shown in the diagram below. This is to be used in real time audio rendering filtering between two FFT filters.

Linear interpolation between each of the FFT's complex coefficients doesn't work as wanted.

FFT plots of two impulse responses
(x = frequency/coefficient, y = frequency magnitude)

+                                   +
|                                   |
|        +                          |
|        |                          |
|        |                          |
|        |                          |
|        |                          |
|        |                          |                   +
|        |                          |                   |
|        |                          |                   |
+--------+------------------+       +-------------------+-------+
           Input A                              Input B

FFT plots of interpolation between A and B (Linear, and Wanted)
(x = frequency/coefficient, y = frequency magnitude)

+                                   +
|                                   |
|                                   |
|                                   |
|                                   |
|                                   |             +
|        +                          |             |
|        |                          |             |
|        |          +               |             |
|        |          |               |             |
+--------+----------+-------+       +-------------+-------------+
Linear Interpolation (A+B/2)            Wanted Interpolation

How do I do the interpolation as described in the Wanted Interpolation graph?


texts: Bilinear transformation of continuous time state space system

texts:
I'm trying to understand the derivation of the bilinear transform for a set of continuous time state-space matrices. I've found plenty of websites which list steps to perform the conversion (here 1 or here 2 or here 3) - but haven't found any derivation or insight into how they came about. When I attempt to do it myself, I keep getting terms in z which I can't seem to get rid of.

To make the problem clear, if I define the bilinear transform as:

$$ s=\frac{\alpha\left(z-1\right)}{z+1} $$

Where alpha is usually given as 2/T.

The transfer function of a continuous-time state-space system can be given as:

$$ H(s) = \mathbf{C}(s\mathbf{I}-\mathbf{A})^{-1}\mathbf{B} + \mathbf{D} $$

How do I derive the expressions for the matrices of the discrete time version bilinear transformed model:

$$ H(z) = \mathbf{C_d}(z\mathbf{I}-\mathbf{A_d})^{-1}\mathbf{B_d} + \mathbf{D_d} $$

When I make the variable change, I keep getting extra terms in z which I don't know what to do with. I'm kinda hoping this is just me failing to have a strong grasp of matrix operations. Any insight (including if I am going about this completely the wrong way) would be really appreciated!

--- Edit with answer.

Thanks to Klaz for posting an answer to this first. I managed to work through the problem today and thought I would post my own alternate working here.

Given the typical state-space equations:

$$ s\mathbf{Q}(s) = \mathbf{A}\mathbf{Q}(s) + \mathbf{B}\mathbf{X}(s) $$
$$ \mathbf{Y}(s) = \mathbf{C}\mathbf{Q}(s) + \mathbf{D}\mathbf{X}(s) $$

Focus on the state update equation first and make the bilinear substitution:

$$ \frac{\alpha\left(z-1\right)}{z+1} \mathbf{Q}(z) = \mathbf{A}\mathbf{Q}(z) + \mathbf{B}\mathbf{X}(z) $$

We manipulate this expression to make the LHS $z\mathbf{Q}(z)$:

$$ \alpha\left(z-1\right) \mathbf{Q}(z) = \mathbf{A}\left(z+1\right)\mathbf{Q}(z) + \mathbf{B}\left(z+1\right)\mathbf{X}(z) $$

$$ z\left(\alpha\mathbf{I}-A\right)\mathbf{Q}(z) = \left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \mathbf{B}\left(z+1\right)\mathbf{X}(z) $$

$$ z\mathbf{Q}(z) = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B}\left(z+1\right)\mathbf{X}(z) $$

This instantly gives us $\mathbf{A_d}$ as:

$$ \mathbf{A_d} = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right) $$

We know that the linear $(z+1)$ term of $\mathbf{X}(z)$ can be applied wherever $\mathbf{Q}(z)$ is used. We apply this in the output expression to give the new output discrete equation as:

$$ z\mathbf{Q}(z) = \left(\alpha\mathbf{I}-A\right)^{-1}\left(\alpha\mathbf{I}+A\right)\mathbf{Q}(z) + \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B}\mathbf{X}(z) $$
$$ \mathbf{Y}(z) = \mathbf{C}\left(z+1\right)\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$
$$ \mathbf{Y}(z) = \mathbf{C}z\mathbf{Q}(z) + \mathbf{C}\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$

We can now pull out $\mathbf{B_d}$ as:

$$ \mathbf{B_d} = \left(\alpha\mathbf{I}-A\right)^{-1}\mathbf{B} $$

Because we know an expression for $z\mathbf{Q}(z)$, we can formulate $\mathbf{Y}(z)$ as:

$$ \mathbf{Y}(z) = \mathbf{C}\left(\mathbf{A_d}\mathbf{Q}(z) + \mathbf{B_d} \mathbf{X}(z) \right) + \mathbf{C}\mathbf{Q}(z) + \mathbf{D}\mathbf{X}(z)$$

$$ \mathbf{Y}(z) = \mathbf{C}\left(\mathbf{A_d}+\mathbf{I}\right)\mathbf{Q}(z) + \left(\mathbf{C}\mathbf{B_d} + \mathbf{D}\right)\mathbf{X}(z) $$

Giving:

$$ \mathbf{C_d} = \mathbf{C}\left(\mathbf{A_d}+\mathbf{I}\right) $$
$$ \mathbf{D_d} = \mathbf{C}\mathbf{B_d} + \mathbf{D} $$


texts: Let's say you have an original image and a version of the same image that may have been convoluted with a Gaussian blur. How could you demonstrate that the Gaussian blur has been applied and calculate the blur radius?


texts:
I have four microphones which I want to use to locate the direction of an audio source. Is there any Python module that implements an algorithm that, given the recordings of the four microphones, will give me a vector pointing to/angles of the audio source?
Thanks in advance


texts: DFT equivalent circular convolution weight matrix with a symmetric filter of length 2K+1

texts:
$\DeclareMathOperator{\diag}{diag}$In a research paper, I read that:

For optimization, the $n\times n$ weight matrix of DFT can be equivalent to circular convolution with a symmetric filter of length $2K+1$. That is the matrix rows are circularly shifted versions of $ h_{-K},...h_{-1},h_0,h_1,...h_K,0,...0)$.

PS: the DFT matrix is
$A = W^{-1}\diag(H_1,...,H_n)W$
, where $W$ is $n\times n$ DFT matrix, and $H_k=\frac{\alpha}{2}+j\frac{\beta}{2}(2\pi f_k)^2$

What exactly does this mean?  How can I get the expression of $h_K$?


texts:
I am sorry if I say any thing silly .Please forgive me.I want to apply a window to input signals.I have two signals each of which contains zeros (462 samples) and non zeros data (38 samples).I want my window to take only the non zeros data (or data above certain threshold).My ultimate goal is to reduce the processing time and resources by taking only non zero data .I will use this data in finding cross correlation and finally the delay between those signals.Can any one comment How I can do that (some hints about matlab code).
Thanks in advance and I am sorry for any inconvenience!


texts:
I'm calculating FFT from microphone input. I notice that lower frequencies always seems to have more power (higher dB) than higher frequencies.

I cut the data into frames of 24576 bytes (4096*6).
Apply Hamming window: input[i] *= (0.54d - 0.46d*(double) Math.Cos((2d*Math.PI*i)/fs));
Run it through FFTW Process1D().
Convert from complex numbers: output[i] = 10.0 * Math.Log10((fout[i * 2] * fout[i * 2]) + (fout[i * 2 + 1] * fout[i * 2 + 1]));
Average out 6 values to get a complete FFT of 4096 bytes.
Paint pretty picture (colors mapped to a palette).

The averaging (pt. 5) is done to lower FFT noise.

As the image shows both with sound and with mic off there is more energy (and more noise) in the lower frequencies. This indicates its more than just a mic/sub problem.

My questions:
 1. Is this expected? Why?
 2. Any standard way of fixing this? Almost looks like some Math.Tan() magic could lift it up where it belongs.

My goal is to be able to identify the top 5 frequencies without lower frequencies winning by default.


texts: How do I convert a filter from Spatial to Fourier domain?

texts:
I have a small 2D Gaussian image filter that I am currently applying as a window by convolution. I would like to do this in the Fourier domain as a single multiplication. I would like to convert my current small filter to the Fourier domain as it is with maximum possible fidelity.

I do not know the size of the image that I shall be applying it on and want to be able to dynamically change it for any image size. How do I go about doing this?

Can I scale my filter to the image size (square to rectangle), FFT the scaled and now distorted and multiply in the Fourier domain?

or

Should I FFT it first, scale up (square to rectangle) and distort the FFT and then multiply it?

I am not convinced that either is a correct approach. Are there any alternative ways?

I am working with this in MATLAB


texts:
Can someone (maybe @hotpaw?) Explain How to use a size N complex DFT to compute a size N/2 real DCT?

Also, I noticed that the DCT matrix (which can be obtained in Matlab by dctmtx(n)) contains dct basis  vectors  that are not whole (integer) frequency numbers e.g. the second basis vector is only 1 half of a cycle (0.5 Hz). Why is that so?


texts:
I appreciate any help: I wanna convert IQ data to its RF data, representing an image. The center frequency and sampling rate are available. What I do is:

First I do up-sampling in MATLAB by applying interpolation along depth. So, the new data has a new sampling frequency higher than the first one.

My questions are: 1) do I need to do low-pass filtering before multiplying the data by the exponential term (that converts IQ to RF)? 2) what is the exponential form that I need to use to multiply the IQ data by?

I was wondering if the length of one of my vectors is, say, 1000, ( Calling it vector A), then how to do zero padding? (I assume a number of zero elements are added to both start and end of this vector and then the Low-pass filtering is done) How to decide about the number of zeros added? Sorry, if these seem to be too simple questions

I've read the documentation I copy here and it is confusing for me: (Page 8 and afterwards)

http://folk.ntnu.no/htorp/Undervisning/TTK10/IQdemodulation.pdf

Thanks


texts:
Sorry if this is not the right forum.
Having a velocity grid in which a wave propagates, and given a certain spatial sampling of the grid, what kind of smoothing/interpolation should be used if the grid is to be resampled (decimated) ?
In the top picture below, I have a velocity grid with a grid resolution of 5 meters. I need to decimate to 25 meters instead. The colorbar shows the velocity variation. The middle picture is the first gradient along the vertical and the bottom picture is the second gradient.
Given waves with central frequency from 10Hz to 30Hz, what kind of smoothing should be applied before decimation to avoid spatial aliasing (least possible) ?

Thanks


texts: Inverse $z$-transform of a transfer function in MATLAB

texts:
I have designed a Butterworth highpass filter (HPF) of 4th order with cutoff frequency high enough to give a gain of $3$ at high frequencies. I want to find the inverse $z$-transform using MATLAB.

L=4;
OSR = 16;
[b,a] = butter(L,0.259,'high');
b = b/b(1);
fvtool(b,a);
H = tf(b,a,[],'variable','z^-1');

iztrans doesn't seem to support tf or zpk.


texts:
While solving Example 4.1 of Signals and Systems by Alan Oppenheim. Example 4.1 is:
$$ x(t)=e^{-at}u(t), a>0$$
and the transform I get is:
$$ X(j\omega)\frac{1}{a+j\omega}, a>0$$

The problem is understanding the sketch:
The sketch is drawn for both the values $+a$ and $-a$> i dont understand why it $-a$ when $a>0$ is defined.


texts: Effect of Gaussian Blur on Different Frequency Components of an Image

texts: Gaussian blur is lowpass filter which means it pass the low frequency components  and stop the high frequency component but the point  I am not getting is that if it is lowpass filter then how it is possible to prevent the edges which is nothing but the "HIGH FREQUENCY COMPONENT"?


texts: Modeling a signal through the environment

texts:
I have a signal that is composed of a sinusoid.  I will send it over the air through radios.  Part of the signal will go through processing on a computer, which means that it will be converted from analog to digital, and back again.

I have two questions.  How can I model any frequency shift that will occur in this system?  And, how can I model how the time of the signal may change?  That is, if the signal occurs in bursts (e.g. 200 ms of tone followed by 100 ms of absence of tone), how can I determine how the time may decrease?  Is the time changed only by the digital part of the system (due to sampling), or might the radio path also affect the duration?

Thank you


texts: Analyzing accelerometer data: appropriate use of rotation matrices?

texts:
I'm working on a project in which a rod is attached at one end to a rotating shaft. So, as the shaft rotates from 0 to ~100 degrees back-and-forth (in the xy plane), so does the rod. I mounted a 3-axis accelerometer at the end of the moving rod, and I measured the distance of the accelerometer from the center of rotation (i.e., the length of the rod) to be about 38 cm.

Below is the plot of both relevant dimensions collected from the accelerometer (after filtering). I am neglecting to include the third up and down dimension as the accelerometer shows a near constant ~1 G, so I think its safe to say its not capturing much rotational motion. Red is what I believe is the centripetal component, and blue is tangential.

Note that I'm also colleccting position and time data, though I don't know if that helps any.

My questions pertain to the following:

1) How do I "combine" these data to fully capture the rational motion taking place? I know it's not as easy as summing the squares and taking the sqrt, but I've been having a really hard time on finding a way to do this throughout all the research I've done.

2) How do I correct for the offset in the data? After some research and asking others, it may have to do with the orientation of my accelerometer (its not perfect). Apparently I should use rotation matrices to transform the measurements to a different coordinate system - one where one axis feels an acceleration of g and the rotation of the rod influences another axis only, but I'm unsure of how to do this.

Can someone help me use rotation matrices on the data? Will using rotation matrices fix the problem listed in number 1 as well? I'd appreciate any kind of assistance I can get. Here are links to both dimensions captured from the accel (the full data):

http://cl.ly/3n3Y2y3j1P19?_ga=1.76624790.2093327149.1426657579

http://cl.ly/3R1d3E1Y2C2x?_ga=1.76624790.2093327149.1426657579


texts:
I am having a very hard time to implement the DCT algorithm. I have quite a few requirements like it has to work with NxN matrix or at least power of 2, it has to be 2D, it has to produce same output as FFTW fftwf_plan_r2r_2d(FFTW_REDFT10) it has to use real data and I need DCT I, II and III. And it has to be fast!

1) I could use FFT to compute DCT and wikipedia mentions this:

"One can also compute DCTs via FFTs combined with O(N) pre- and post-processing steps. In general, O(N log N) methods to compute DCTs are known as fast cosine transform (FCT) algorithms."

What are those 2 steps?

2) If I use Apples vDSP library is that even a good idea? since it needs an array twice as big as the DCT array (2N with radix2). And also vDSP is 1D only so I would end up doing it for each row and column to get a 2D result.

3) Is it possible to use an algorithm that works on 8x8 blocks and adapt it for my needs?

I am very new to learning about DSP so any help is appreciated! Thanks!

ps: Does anyone have any sample code for what I need? would help me a lot


texts: Contradiction while using the convolution sum for a non-LTI system

texts:
In a recent quiz, we were given the following problem:

  The cascaded LTI systems $\mathcal{T}_1$ and $\mathcal{T}_2$ respectively have impulse responses $h_1[n]=\delta[n+3]$ and $h_2[n]=\delta[5−n]$. What is the output when the input is $x[n]=n$, i.e., find $y[n]=\mathcal{T}_2\{\mathcal{T}_1\{x[n]\}\}$.

I am primarily confused about the "LTI-ness" of the system and at a contradiction; given the second stage to be LTI (assuming that for a system to be LTI as a whole, all its stages/subsystems must also be LTI). Thus, we can say that a sub-system with  $h \lbrack n \rbrack = \delta \lbrack 5 - n\rbrack$, i.e. $y \lbrack n \rbrack = x \lbrack 5 - n\rbrack$ should be LTI, which already seems to be false. To prove that, I considered the following convolution sum: $$\tilde{y} \lbrack n \rbrack = \displaystyle{\sum_{k = -\infty}^{\infty} x \lbrack k \rbrack \delta \lbrack 5 - n + k \rbrack} = x \lbrack n-5 \rbrack \ne x \lbrack 5-n \rbrack = y \lbrack n \rbrack$$

And this is a contradiction since, after the convolution, the output doesn't match the original output with which we started.

Also, I realized that interestingly, the convolution always gives a result that corresponds to an LTI system (as in the case above too: $y[n]=x[n−5]$ is LTI whereas $y[n] = x[5-n]$ is not).

Also, since $\delta[n]$ is even, any non-LTI system's impulse response, for example, $\delta[1−n]$, will equal $\delta[n−1]$ which corresponds the impulse response of an LTI system. This explains why I am getting LTI characteristics after the convolution. This is interesting too since in either way the impulse response implies delaying the signal by 1 (in the current example).

And so did the impulse response $\delta[5-n]$ stated in the quiz question just qualitatively imply a delay of 5, and is, technically, not the exact description of the underlying system?

So, in summary, I have the following two doubts:

Is the convolution sum only true for the output relation of an LTI system? If so, then can it be used to prove non-LTI-ness of a system in the same manner as above?
Is the quiz question incorrect?

Any help will be greatly appreciated!


texts: How to shift the center frequency of an IIR filter

texts: I have a pre-designed in-series 14th-order IIR filter in second-order sections. Problem is, when I simulate and measure the behavior in my codec, there is a frequency shift downward (as compared to the MATLAB prediction). Is there a simple way to alter the coefficients in such a way to shift the center frequencies of each band? Thank you.


texts:
I was just learning about psychoacoustic audio data compression and came across some sources suggesting that a strong tonal component (masker) that occurs within a certain short timespan after a weaker component (maskee) has the effect on the brain that, it will be unable to perceive the maskee.

What I find really hard to digest is how can it affect the causality of time? Is it because the maskee signal gets picked up by the ear but gets dropped before it gets perceived in the brain due to the masker signal happening before the communication/processing delay between the ear and the brain?

Any thoughts on this?


texts: How to derive the transmitted signal from a constellation diagram

texts:
Let us assume we have $16$ symbols to transmit.  We can represent these $16$ symbols by $16$ unique signals and transmit.  If at the receiver we can identify these $16$ signals correctly, we have identified the symbols transmitted correctly.  I have used qammod([0:15],16,0) in MATLAB and got the output that is attached here.

>> a = qammod([0:15],16,0); a'

ans =

  -3.0000 - 3.0000i
  -3.0000 - 1.0000i
  -3.0000 + 1.0000i
  -3.0000 + 3.0000i
  -1.0000 - 3.0000i
  -1.0000 - 1.0000i
  -1.0000 + 1.0000i
  -1.0000 + 3.0000i
   1.0000 - 3.0000i
   1.0000 - 1.0000i
   1.0000 + 1.0000i
   1.0000 + 3.0000i
   3.0000 - 3.0000i
   3.0000 - 1.0000i
   3.0000 + 1.0000i
   3.0000 + 3.0000i

$0,1,\ldots,15$ represent the $16$ symbols. $16$ represents the number of symbols or the number of unique signals - sinusoids with different amplitudes and phases - to be transmitted.  $0$ represents the offset phase.

Now I want to interpret the output.  I actually got $16$ complex numbers.

What do they represent?
What are the magnitude, phase, real part and imaginary part of these complex numbers that they represent?  I guess they represent the parameters of sinusoidal signals.
In both the real and imaginary parts we see $-1, -3 , 1$ and $3$.  We do not see $-2, 0$  and $2$.  Why?


texts: I have a collection of data: number of hours per day, over the course of several years.  I want to create an FFT of this data.  Since my sampling rate isn't seconds, but days, I think I have two choices when constructing the FFT.  One, I can specify that the sample rate Fs=1/86400 samples per second (86400 seconds in a day.)  But I think that's not right.  So maybe I can: two, say Fs=7 samples per week, then my resulting FFT has a scale, instead of cycles-per-second, of cycles-per-week.  Is this reasonable, or is there a better way to go about doing FFT for daily-periodic data?


texts: Decorrelating Stationary Colored Gaussian Noise -- Effect On The Desired Signal

texts:
So given stationary colored gaussian noise $\mathbf{n}$, I know that I can decorrelate it by first finding it's autocorrelation $R_{nn}$ and performing $R^{-\frac{1}{2}}_{nn} \mathbf{n}$.

In practice of course, I need to estimate $R_{nn}$, which I can do via averaging $\mathbf{r}_{nn}$ or by fitting it to an AR model.

So, suppose I have a good estimate $\hat{R}_{nn}$, for a finite data record of length $N$. Such that the resulting spectrum of the noise frame is white enough, according to some metric.

Now, suppose I have a signal $\mathbf{x} = \mathbf{s} + \mathbf{n}$. I want to decorrelate the noise, such that after the decorrelation, I get a signal in white noise. From what I have read, the way to do this is to also apply $R^{-\frac{1}{2}}_{nn} \mathbf{x}$.

However, won't this distort the desired signal $\mathbf{s}$? Can anyone point me to some references regarding this?

Also: How else can I decorrelate colored noise in order to get my signal embedded in white noise?


texts: Autocorrelation function of a discrete signal

texts:
For a certain experiment of which I have data, I want to calculate the autocorrelation of the measurements. For this I used the expression from Wikipedia, which reads

Here, $n$ is the dataset's size, $0 < k < n$, $\sigma$ is the standard deviation and $\mu$ is the mean.

Now, I've read the descriptions on the page, but I don't fully understand all of the terms in there. I understand why we divide by $(n-k$), because the number of datapoints used constantly changes (you have a finite signal, after all). However, I don't see why we divide by the variance. Moreover, I also don't know why we subtract the mean. Doesn't this lead to problems, if for example for n = 50 and k = 49, we want to multiply $(x_1-\mu)(x_{50}-\mu)$, and $x_1 > \mu$ but $x_{50} < \mu$? I don't really see what a negative correlation would entail.

For my data, which is only 0's and 1's, this is even more problematic. Using the above formula the correlations quickly decay to zero, but at the end they violently start to spike, because the contribution from the mean becomes so big, relatively. In my data the mean is 0.5 and the variance is 0.25, so for example, my final point $R(49)$ = $\frac{(0-0.5)(1-0.5)}{0.25}$ = -1

But -1 is huge, especially when 'previously' the correlations had died out to 0. How does one deal with this?

An illustration of this is


texts: Why normalize the data set before applying Direct Linear Transform

texts:
Direct Linear Transform (DLT for short) is a method of homography estimation, it solves the overdetermined linear system via SVD $$Ah=b$$to find a solution $h$ under constraint $\|h\|=1$. Actually it finds the least square solution which minimize $\|Ah - b\|$.

I understand the basic idea of this algorithm, but it is recommended to normalize the data set before applying DLT on it, and here is a intro about how to do the normalization. It is lectured that data normalization is important to DLT, without normalization the results from DLT is not stable.

I wonder why? Just because DLT involves solving the linear system using SVD and $A$ might be singular?


texts: How to Apply Least Squares Estimation for Sparse Coefficient Estimation?

texts:
The model is expressed as, $$y(n) = \sum_{i=0}^{p-1} r(i) x(n-i) + v(n) \tag{1}$$ where  $\mathbf{r} = [r_1,r_2,\ldots,r_p]^T$ is the sparse channel coefficients of length $p$, $\mathbf{x} = [x_1,x_2,...,x_{n-p}]^T$ is the one dimensional input and v is additive White Gaussian noise of zero mean. y = filter(.) command is used to model the above equation and thus creating an FIR filter or a moving average (MA) model. The order of the MA model is p=3.

So, $y = [y(1),y(2),....,y(100)]$ is a vector of 100 elements.  I am generating noise of variance 0.1 I want to estimate the sparse channel coefficients using LASSO. As there are p channel coefficients, I should get p estimates.

According to the equation of LASSO, ||rx - y||_2^2 + lambda * ||r||_1 I am estimating the sparse coefficients, r. As the true coefficient array contains p elements, I should get p estimated elements. I am not quite sure if this is the way to do. I have not found any example on LASSO being applied to univariate time series model such as ARMA. I don't know how to estimate the sparse coefficients using the appropriate algorithm and need help.

The first part of the Equation : $||rx - y||_2^2$ is  a least squares formulation which I can solve using Least Squares Approach. In order to implement LS, I have to arrange the input in terms of regressors. However, if the coefficients, $\mathbf{r}$ are sparse then I should use LASSO approach. I have tried using Matlab's LASSO function. For LASSO, I rearranged the input data $x$ in terms of regressors, but I don't know if this the correct approach.

I need help. Is there an approach to include the sparsity term in the LS?

Please find below the code for LASSO using Matlab function. As a toy example I am just assuming model order to be of lag 3 but I know that LASSO can be applied efficiently to a large model. I can test for larger order MA model having lag > 3.

% Code for LASSO estimation technique for
%MA system, L = 3 is the order,

%Generate input
 x = -5:.1:5;

r = [1    0.0   0.0];% L elements of the channel coefficients
%Data preparation into regressors
X1 = [ ones(length(x),1) x' x']; %first column treated as all ones since    x_1=1

y = filter(r,1,x); % Generate the MA model
[r_hat_lasso, FitInfo] = lasso(X1, y, 'alpha', 1, 'Lambda', 1, 'Standardize', 1);

OUTPUT :

The estimates returned are r_hat_lasso = 0, 0.657002829714982, 0

Question : This differs very much from the actual r.
Is this wrong?


texts: Differences between low-pass, band-pass, notch filters

texts:
I'm studying image processing and I know that there are different filters you can apply. I don't understand the difference between the low-pass/bandpass/notch filter (or high-pass/band-reject/notch). I tried to draw them but I don't know if I did it correctly:

From what I have realized it is that the difference between a low-pass filter and a band-pass is that the first is symmetrical with respect to the y axis while the second is not. Therefore, a low-pass filter passes the lower frequencies and resets those higher than a certain threshold (cutoff frequency). Instead, a band-pass filter passes the frequencies included in a certain range and cancels the others. It's correct?

But I have not really understood the purpose of the notch filters. I realized that there are notch band-pass and band-reject notch. In addition the book I'm studying (Digital Image Processing - Rafael C. Gonzalez) says that a notch filter is composed of the product of two high-pass filters. Just I don't understand how they work, how they are made and what they do.

The book also talks then about optimal notch filters that is used where there are several interference components in the image in the Fourier domain. Again I haven't understood how it works.

Could someone help me? Thank you very much!

So the difference between a notch filter and a band-pass is that I use the notch when I don't know what are the frequencies to pass (or delete) but I get it visually analyzing the image of the Fourier transform, right?

However, this is what there are written in my book:


texts: Statistics of the output of a high bitrate signal filtered through a narrow band low pass filter

texts: I was wondering what would be the statistics of the output of the filtering of an extremely wideband signal (e.g. 10 Gbps random data stream) through a narrowband low-pass filter (e.g. 10 MHz LP). I have the feeling that the statistics of the output would be somehow Gaussian due to the central limit theorem, but I cannot give a clear explanation to this.


texts:
I prepared a Matlab script, which re-samples spectrum from 24kHz to 20kHz. When I used up-sampling, a fir filter, and then up-sampling again at the end, I get an additional peak at 4kHz which is not expected. When I checked via the resample function from Matlab paltette, there is no additional peak. What did I do wrong? What do I have to add to my script to avoid this kind of situation. Thank you in advance

clc;
close all;
clear variables;
Fs1 = 24000;
TimeStep1 = 1/Fs1;
Lenght = 10000;
time1 = (0:Lenght-1)*TimeStep1;
InputSignal = sin(2*pi*8000*time1);
NFFT1 = 2^nextpow2(Lenght);
Y1 = fft(InputSignal,NFFT1)/Lenght;
f1 = Fs1/2*linspace(0,1,NFFT1/2+1);
figure
plot(f1,2*abs(Y1(1:NFFT1/2+1)))
title('Original Spectrum')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
Fs2 = 20000;
T2 = 1/Fs2;
nGcd = gcd(Fs1, Fs2);
upY = upsample(InputSignal,Fs2/nGcd);
firCoeff = fir1(255, 0.9, 'low');
filteredY = filter(firCoeff,1,upY);
dataOut = resample(InputSignal,Fs2,ceil(Fs1));
L2 = length(dataOut);
t2 = (0:L2-1)*T2;
N=2;
lenght = (L2/N);
NFFT = 2^nextpow2(lenght);
Y = fft(dataOut,NFFT)/lenght;
f = Fs2/2*linspace(0,1,NFFT/2+1);
figure
plot(f,2*abs(Y(1:NFFT/2+1)))
title('Resampled Spectrum')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')
figure
dataOut2 = downsample(filteredY,Fs1/nGcd);
Y2 = fft(dataOut2,NFFT)/lenght;
plot(f,2*abs(Y2(1:NFFT/2+1)))
title('UpFIRDwn Spectrum')
xlabel('Frequency (Hz)')
ylabel('|Y(f)|')


texts: Is it correct to assume that $(-1)^n = \cos(\pi n)$ while computing the DTFT of $(-1)^n$?

texts: $\cos(\pi n)$ fluctuates between $-1$ and $1$ depending on the values of $n$, and it would be the same as computing it with an exponential but the problem is that I just get part of the right answer...


texts: Finding the memory polynomial model coefficients

texts:
I want to build a memory polynomial model given that I know the corresponding output signals for known input signals.
What is the easiest way to find the coefficients of it?
And how do we take the memory and order into consideration while we build it.

This is what I know. if A and B are input and output signals of length n, then A_inverse x B will give me n by n  matrix.
But how can I introduce memory and order?, if i add any extra columns to A_inverse x B, it will change the dimensions to n by m where m > n. and i can no more operate it on A.
 'x' means matrix multiplication.


texts: Compensating for group delay in FIR with symmetrical coefficients

texts:
Suppose I have a sinusoidal signal of 50Hz sampled at 6.4kHz and I wished to filter it at:

Band pass: 160Hz, Bands stop: 800Hz, Apass:0.1db, Astop: 106db

For a FIR filter order of 40 I would get 41 number of taps. The group delay would thus be (41-1)/2 = 20.

Am i wrong in thinking that if I extended (/zero padded at the end of) the data block size to be filtered by the estimated group delay above, I would be compensating for the group delay?

For example if I had 128 samples (corresponds to a full sine cycle at 50Hz-Fs=6.4Khz) to be filtered and added the group delay 128+20=148. Now, if I were to filter these 148 samples with the above stated filter, would I not get a signal with no group delay?

Thanks
Alex


texts: python3: speech recognition - calculating features of a WAV file

texts:
I am supposed to do this as a school project, but I am kind of lost when it comes to signals and would appreciate your help. Also, English is not my first language, and I will be translating the assignment mostly literally, as I am not familiar with a lot of terms used in signal processing, so pardon me for any mistranslations, please.

I have WAV files with these characteristics:

mono
16000 Hz
16-bit precision
16-bit Signed Integer PCM encoding

Right now I am supposed to calculate a vector of features for each frame (sample?) of each WAV file.

They even suggested a solution. They recommended to do this using a linear filter bank, which produces 16 coefficients by summing coefficients of a logarithmic power spectrum. There's an exact formula, but it's irrelevant for this question, I think.

What can I do to get an array with values of this logarithmic power spectrum (is that the right translation)?

I have imported these libraries and can freely use any of its functions: IPython, numpy, matplotlib, soundfile, scipy.signal

Thank you


texts: Converting Gamut3D from RGB colorspace to LAB colorspace

texts:
I have this code for Gamut3D in RGB, but I'm looking to convert it into LAB colorspace but I'm not really sure what changes I should be making.

rgbImage = imread('pill1.jpg');
% Get the dimensions of the image.  numberOfColorBands should be = 3.
[rows columns numberOfColorBands] = size(rgbImage);
% Display the original color image.
imshow(rgbImage);
title('Original Color Image');
% Enlarge figure to full screen.
set(gcf, 'units','normalized','outerposition',[0 0 1 1]);
% Extract the individual red, green, and blue color channels.
redChannel = rgbImage(:, :, 1);
greenChannel = rgbImage(:, :, 2);
blueChannel = rgbImage(:, :, 3);
% Construct the 3D color gamut.
gamut3D = zeros(256,256,256);
for column = 1: columns
        for row = 1 : rows
                rIndex = redChannel(row, column) + 1;
                gIndex = greenChannel(row, column) + 1;
                bIndex = blueChannel(row, column) + 1;
                gamut3D(rIndex, gIndex, bIndex) = gamut3D(rIndex, gIndex, bIndex) + 1;
        end
end
% Get a list of non-zero colors so we can put it into scatter3()
% so that we can visualize the colors that are present.
r = zeros(256, 1);
g = zeros(256, 1);
b = zeros(256, 1);
nonZeroPixel = 1;
for red = 1 : 256
        for green = 1: 256
                for blue = 1: 256
                        if (gamut3D(red, green, blue) > 1)
                                % Record the RGB position of the color.
                                r(nonZeroPixel) = red;
                                g(nonZeroPixel) = green;
                                b(nonZeroPixel) = blue;
                                nonZeroPixel = nonZeroPixel + 1;
                        end
                end
        end
end
figure;
scatter3(r, g, b, 3);
xlabel('R' );
ylabel('G' );
zlabel('B' );


texts: I've been studying communications and i see that the book used the symbol $u_{-1}(t)$ for unit step function , why ?


texts: What unit is the z axis of a power spectrogram of audio?

texts:
the question might seem dumb but since I tried it out and got values well above 2000 (watts?) it got me thinking... Check this out:

winLen = 10*1e-3;
fftOverlap = 0.5;
[x, Fs] = audioread('speech_5.wav');
x = x(:,1);
hopSize = winLen * Fs;
nfft = hopSize / fftOverlap;
overlap = nfft - hopSize;
w = hann(nfft);
[S, f, t] = spectrogram(x, w, overlap, nfft, Fs);
S = (abs(S)).^2 + eps; %make it a power spectrogram

If you look into S (the spectrogram) now, you see a lot of values above 1000, in my case (this is clean speech) even 6000. It got me thinking, what kind of unit this is? Intuition says its watts, but 6000 watts of power in a certain frequency in a 10ms time frame seems very unlikely to me... I mean my hair dryer has 2000 watts... Can someone please help me out here?

Additionally, if you see anything in my code that could be improved, please share! I will be more than thankful for any input!


texts: What features of using Fourier transform method in design FIR filter

texts: I want to know when the  Fourier  transform method is best for design  FIR filter


texts: Coherent Sampling And The Distribution Of Quantization Noise

texts:
I have a question concerning the distribution of noise in the frequency domain in case of Coherent Sampling. I read up about Coherent Sampling and understood that in order for a frequency $f_{in}$ to be located on one single bin of the DFT, the period or a multiple of it must fit exactly into the data window of the DFT, meaning that the following must be true:

$$
\frac{f_{in}}{f_s}=\frac{N_c}{N}
$$
$N$ ... Record Length, $N_c$ ... Number of cycles in data window

Furthermore in the slides of a course about Mixed-Signal Processing I am attending it is mentioned that $N_c$ and $N$ need to be coprime in order for the quantization noise to be distributed uniformaly. Otherwise the quantization noise would show a discrete spectrum. Why is that the case?


texts: IIR biquad coefficients for negative frequencies in IQ data?

texts:
Assume I have a (single real float channel) biquad IIR filter implemented in C, and some coefficients for a bandpass filter (right out of RBJ's cookbook).

If I convert the input, output, data path, and internal arithmetic of this IIR filter to use complex data types and multiplies, can I use the resulting code to filter a complex (IQ) data stream?  If so, how might I modify the biquad IIR coefficients to filter just the negative frequency versus the positive frequency band in a complex (IQ) data stream?


texts: If I define a channel to have complex Gaussian coefficients, do I still need to add noise to the signal passing through this channel using the command "awgn" or would that be redundant?


texts: How to determine windows size of moving average filter

texts:
I've recorded acceleration signal data using an accelerometer attached to vehicle , and I want to calculate the displacement of signal data using double integration or any other methods. I read many studies that perform moving average filter to acceleration data to reduce the amount of error caused by the suspension system of the vehicle.
How I can determine the ideal windows size of moving average filter for this purpose?


texts: How do I compare two resultant images from different binarization algorithms?

texts:
I have a binarized image to which i add noise and then try to filter out the noise using various thresholding algorithms such as otsu and niblack.How can i compare the resultant image with the original image so as to find the percentage error that exists between the two??

the original image is as such:

and the resultant image is:

I need a way to find the percentage error that is present.


texts: Spectral window when FFT is not a power of 2

texts:
I am trying to better understand zero padding and the FFT. To do so, I am attempting to plot the spectrum of the Hann window function, i.e.

$$ w(n) = \frac{1}{2}\left(1-\cos\left(\frac{2\pi (n+1/2)}{N}\right)\right) $$

To accomplish this, I did the following: 1) I created a series of points $n=0,1,...,49$ that I substituted into the hamming window function to get $w(n)$, 2) I zero-padded the data to $M=512$ to get $x(n)$, 3) I used the FFT function in R to compute the Fourier transform of $x(n)$, 4) I calculated the magnitude squared of the Fourier transform obtained in part 3). This gave me the following result:

I then repeated the above steps, but with $M=500$ instead of $M=512$ (i.e. I zero-padded to 500 instead of 512). This gave me the following result:

Why are these two spectra so different? Everything that I have read about the FFT seems to say that I should choose $M$ to be a power of 2 in order to reduce computation time. But I haven't found anything that explains why the spectrum should change as it does.

Any help would be appreciated... thanks!


texts: Delete frequency components with FFT and IFFT

texts: I have a signal in Matlab defined by signal and t. As it is a noisy signal, I want to delete every component below a certain given frequency "x" Hz with the FFT. How could I do it? I guess the procedure is FFT the signal, delete those components and then apply the IFFT, but I don't know how to implement it. Any help is appreciated!


texts:
When applying Gaussian filters close to the borders of an image, algorithms usually zero-pad or replicate/mirror/wrap the borders. This is not good enough for my case, so I wonder if there is something better out there.

As the idea of smoothing is to absorb the information of surrounding points, I wonder if it would make sense to use a different kernel in the border cases. This kernel would only take values from available surrounding points, e.g. if my kernel is something like:

\begin{bmatrix}
    0.0113 & 0.0838 & 0.0113 \\
    0.0838 & 0.6193 & 0.0838 \\
    0.0113 & 0.0838 & 0.0113
\end{bmatrix}

To calculate a pixel on the top margin (but not close to the corner), I would use

\begin{bmatrix}
    0.0938 & 0.6933 & 0.0938 \\
    0.0126 & 0.0938 & 0.0126
\end{bmatrix}

The latter matrix was calculated by cropping out the part of the kernel that would be out of the image (the top row) and renormalizing. Is there something wrong with this approach, or an alternative?

Here's some code to explain my idea better:

kernel = [0.0113 0.0838 0.0113; 0.0838 0.6193 0.0838; 0.0113 0.0838 0.0113]
kernel2 = [0.0938 0.6933 0.0938; 0.0126 0.0938 0.0126];
M = [1 2 3; 4 5 6]; % my image

% application of kernel to point 2 in M would provide the value:
res1 = sum(sum( kernel(2:3,:).*M )); % 0 padding, equals 2.1058
res2 = sum(sum( kernel.*[1 2 3; M] )); % row replication, equals 2.3186

% application of kernel2 to the same point:
res3 = sum(sum( kernel2 .* M )); % equals 2.3568

Note that in using kernel2 I don't extrapolate a row of [0 0 0] or [1 2 3], which would underestimate the point.


texts: Time length difference estimation between a waveform and its stretched and a bit distorted copy

texts:
I'm trying to solve the following problem and will appreciate your help.
In the example I have two signals. The 1st is the original signal and the 2nd is a some slightly stretched version of it. The 1st signal I segment into N overlapping windows and for each window I want to find its new length in the 2nd signal in order to know eventually the difference.

As I see, using cross-correlation is not suitable here, because the segmented waveform in the 2nd signal is not the same length as the original.

Using peak detection works only for part of the signal since close to the end the peaks become more and more distorted.

Will be happy for any help.


texts: How to achieve Shannon's capacity if coding introduces a spectral efficiency loss?

texts:
I am a little bit confused concerning different Shannon's notions :

Over an AWGN channel Shannon's capacity is in my knowledge : $R_b<C=B.\log(1+SNR)$.
Is Shannon's "limit" a different notion than that ?

How to convert this equation to a $BER=f(E_b/N_0)$ chart ? for example Here there is a Shannon limit plotted in fig.2, how is it done ?

From the equation stated above one can derive the maximum spectral efficiency : $\eta = R_b/B <\eta_{max} = C/B = \log(1+SNR)$ is that correct ?

If a channel code introduces redundacy bits it introduces a inevitably a loss in $\eta$, so what does Shannon mean when he states that to achieve capacity one has to tend code length to infinity ?  one of my proffessors once said that the LDPC is an example of such a code because it is very long and thus achieves capacity (almost), where is spectral efficiency loss gone ?


texts: How to find out amplitude of harmonic after diff type of FT?

texts:
I am new in DSP, so I apologies if it's dump questions but I didn't find needed information on the Internet.
So the questions is:
Assume that we have simple harmonic signal described by next function $$x(t) = A \cos(2 \pi f_0 t)$$ After Fourier transformation amplitude of initial harmonic in spectrum changed as next:

For FT: $A/2$
For DFT and FFT: $(A \times N)/2$ where $N$ is the number of samples

Am I right about the last one for discrete transformations? And what about DTFT what is amplitude at frequency $f_0$ equal to? Is it equal to $A/(2T)$ where $T$ is sampling period of the signal.


texts: Characterization of transfer functions with no local peaks

texts:
Assume that you are giving an arbitrary amplitude frequency response $A(\omega)=|H(j\omega)|$

Is there a characterization that ensures that $A$ is monotone? i.e, $A$ has a global maximum at $\omega=0$ and does not have a local peak at some other frequency?

thanks


texts:
Does anybody know why we need 2 different tone signals: one high frequency and one low frequency to create a tone burst in Dual-Tone Multi Frequency signalling?

For example, to create a tone burst of 3, we need a high frequency $1477 \mathtt{Hz}$ and a low frequency $697 \mathtt{Hz}$. Does anyone know the reason for that?


texts: The size of an FIR filter for PDM-PCM conversion

texts:
I have been trying to understand how the conversion from PDM to PCM happens so that I could interface with a microphone.  I have been doing some research online and I found that the most common way to do this is to use an FIR filter and a decimator.  I am having some trouble trying to determine the size of the FIR filter that I should use.  Since in this case the FIR filter is using  binary numbers, the output is affected by the 1’s density.  If the coefficient of each tap is the same, then the number of possible values of each average is equal to the number of taps.  If for example we have a 4 tap filter then 0101  would have the same results as for example 0110 since both values have the same density of 1’s.  In this case the could be 5 different densities of 1 (0 – 1, 1-1, 2-1, 3-1, and 4-1) which can be stored in a 3 bit word.  Using the same logic an FIR filter with 65534 taps would be needed to produce a 16 bit result. As already mentioned, in this case I am assuming that the coefficient are equal, example 1/N.  I understand that the if the coefficients were not the same, then the results would be different because then the order of each bit would also make a difference but in this case, the order of the bits is irrelevant therefore a 1010 would yield the same results as 0110.
I hope that I have explained the problem that I am facing well.  Can somebody please explain how a 16 bit resolution is obtained using FIR filter because to build a 65534 tap filter it becomes impractical.


texts:
I am trying to extract only urban aeas from satellite imagery Using Matlab. I have applied some edge detectors like canny and also used morphological functions and I am able to extract all the boundaries ,but not able to seperate them from various urban aeas like roads, trees and waterbodies.

Please help i can share the code if any one wants


texts:
If we subtract the mean of image from the DCT coefficients of an image

what does it signifies?
and what type of distribution does it follow?

Example in MATLAB:

i=imread('cameraman.tif');
c=mean2(i);
d=dct2(i);
f=c-d;

what does f signifies?


texts: Why integrate over $2\pi$ in inverse DTFT?

texts:
In DTFT of a signal, the spectrum of a sequence is periodic with period $2\pi$ and all the information needed for derivation of the original signal from its spectrum is contained in $\pi <\omega <\pi$.

But , why do they integrate over $2\pi$ and not from $-\infty$ to $+\infty$ , as in continuous Fourier transform in inverse DTFT? After all, when you derive a transform of a signal, the inverse of it is not arbitrary and to recover the original signal, you can't integrate over an other  interval because all the information is contained in there.


texts:
Given the absence of information about the algorithms of spectral whitening functions in software packages:

How many ways do you know for doing spectral whitening?
What are the differences?


texts: Slip rate for sliding correlator channel sounder in GNU radio

texts:
I have a spread spectrum channel sounder running the following transmitter and receiver:
[![tx][1]][1]
[![enter image description here][2]][2]
Imm not sure how to modify the slip rate between the generated pn sequences at the transmitter and receiver. Is it just tied to sps specified at the usrp?


texts: How does SC-FDE modifies $E_b$ (in AWGN function)

texts:
In lot of matlab scripts of SC-FDE simulations I see $E_b = N$ (being $N$ the length of the block to be transmitted through the channel and Eb energy per bit) to generate the AWGN. I would like to understand why is that, and more specifically:
Assuming, instead of a block length $N$, I have 2 sequences length $N/2$. I FFT() each one, concatenate the frequency-domain sequences and send the result (of length $N$) through the channel. What would my $E_b$ in AWGN function be?


texts: MRI reconstruction using windowing based apodization

texts:
I am trying to apply windowing to do MRI reconstruction. I have a 256 point one dimensional k-space samples which look as follows:

Here, the real part(green) is flat around 0 and the blue curve is the imaginary part. Now, what I would like to do is perform some apodization. For this, I construct a hanning window as follows:

from scipy import signal

win = signal.hann(256)
# sig is the
sig = signal.convolve(sig, win, mode='same') / sum(win)

The hanning window is the same size as the k-space vector (256). Next, I convolved the k-space signal with the hanning window, which looks like:

Next, I performed the inverse fft to reconstruct the image (which should have near uniform intensity everywhere) using:

from scipy import fftpack

rec = np.abs(fftpack.ifftshift(fftpack.ifftn(sig)))

This results in the recon image as:

I feel I have done something wrong as I was expecting to have a fairly uniform image rather than a spike. Do I need to do something post reconstruction to account for the fact that I have performed some windowing before the reconstruction?


texts:
I'm trying to figure out which of the resampling methods is best in my job. Comparing the original image vs. bilinear, bicubic, bicubic spline, lanczos. For now I have considered metrics like MSE RME SSIM and UQI, but I read in some articles that the SRE (signal to reconstruction error ratio) is also valid for evaluation. Should I consider it?
If yes, how is it implemented on python?

Finally, the mean square error seems too high.

                                MSE            RMSE      PSNR       SSIM         UQI

NORMAL vs BILINEAR        1842447.8454      1357.3680   21.4111     0.04        0.8548

NORMAL vs CUBIC           1853309.4604      1361.3631   21.3855     0.05        0.8541

NORMAL vs CUBICA SPLINE   1834833.0737      1354.5601   21.4290     0.03        0.8552

NORMAL vs LANCZOS         1855455.1160      1362.1509   21.3892     0.05        0.8540


texts: Deconvolution of a 1D Time Domain Wave Signal Convolved with Series of Rect Signals

texts:
I have a synthesized signal (the bottom of the following figure), which is the convolution of the input signal (at the top) and the objective function (in the middle). The intention is to retrieve the objective function from the convoluted signal, when the input signal is known. From practical point of view, it seems to be an ill-conditioned problem, but I'm curious to know the expert opinions on how far one can get, and with what SP toolset. Specifically, I would like to know: 1). comparing the first wave packet and second, how to tell whether the first square in the objective function is one with shorter time span but the same amplitude, or with the same time span but smaller-amp; 2). for the 4th square, how to tell it's not two separate short squares rather than a single long one. I've also attached the Matlab codes I used to generate the figures. Thank you very much.

clear; close all;

%% 10MHz incident signal;
fs = 1000e6;
f = 10e6;
wavelength = fs/f;
sig = wavemaker(3.5, f, fs);
figure; subplot(3,1,1);
plot(sig, 'LineWidth', 1); title(strcat('input signal, wavelength =',num2str(wavelength),' data pts'));
axis([0,3300,-2,2])

%% Sparse signal;
N = 3000;                    % N : length of signal
s = zeros(N,1);
k = [50:(50+wavelength*0.1) 500:(500+wavelength*0.6) 1200:(1200+wavelength*1.6) 2200:(2200+wavelength*4.6)];
s(k) = 1;
subplot(3,1,2);
plot(s, 'LineWidth', 1); title('distribution: objective function');
axis([0,3300,0,2])

%% convoluted signal
y = conv(sig,s);
subplot(3,1,3);plot(y, 'LineWidth', 1);hold on;
plot(abs(hilbert(y)), 'LineWidth', 1);
title('convoluted signal between input and districution');
xlim([0,3300])

function x = wavemaker(nCycles, fc, fs)
% function to generate wave packet;
nSample = round(fs / fc * nCycles);
ts      = 1 / fs;
T       = ts * nSample;
t_max   = ts * (nSample-1);
t       = 0: ts: t_max;

x       = sin( 2 * pi * fc .* t);

x = x.*hanning(nSample)';
end


texts:
I designed the filter below on paper:

FIR
low pass
type 1
length 3
alpha = 6
beta = -12
allow freq. <= 0.1rad/sample
blocked freq. >= 0.4rad/sample

To confirm that my calculations are correct I was planning to design the filter on MATLAB. Is there a function to design this filter, apply an input frequency and check the output result?


texts: FM demodulation amplifier with saturation

texts:
I'm trying to understand how FM demodulation works. I read a lot on the internet and as far as I understand, the main process is to extract the phase difference (that is also the output) between a quadrature signal and the input signal, and that use this difference to keep those two signals locked in phase.

However the output signal amplitude depends a lot on the input signal strenght. As I understand, the input of the phase locked loop should be constant in amplitude.

Besides on this article on wikipedia it's written that FM recivers do not need a AGC, but it's enough to use a "high-gain IF amplifier which is intentionally driven into saturation". How can this be right? I mean.. If the amplifier is driven into saturation, the input wave will look more like a square wave and all the quadrature equations won't be true anymore! Can anybody explain me this step?

EDIT 1

I think that the input should be constant in amplitute because(I'm thinking about a PLL):

$ FM_in(t) = A sin( W_c * t + \phi(t) ) \\
 FM_in(t) * cos(W_c*t) = {A \over 2}(sen(\phi(t)) +  cos(2W_c + \phi(t)) ) $

So filtering out the high frequency cosine, I get:

$ DEM_Out(t) = {A \over 2} sen(\phi(t)) \simeq {A \over 2}\phi(t) $

And so far so good, but to keep in lock the loop I should feed this value to the VCO, and if A doubles(as I understand), the phase shift of the VCO doubles, too so that the loop won't be in lock anymore. I hope I have explained what I mean!


texts: Frequency response of a rolling linear regression

texts:
I am looking for a way to characterise the frequency response of the slope from a linear regression. We are exploring the effect of window length of the regression to the magnitude of the slope of the regression and find that for our data (discrete time series data), the slope is positively related to the window length.

The slope can be thought of as the expected change of our output variable for a given input variable over a time period equal to the window length. We find that for a short window length, this slope is very low whereas for a long window length the slope is high. We interpret this to mean that the system responds to changes in the input in the long term but not in the short term. To me, that sounds a lot like a low pass filter killing off the fast changes in the input but passing the slower (or rather, more sustained) changes.

The slope, or beta, of the regression for a rolling window is, in a way, similar to a moving average and it is possible to characterise the frequency response of a moving average filter as a function of its window length as shown here. I was wondering if anyone knew of a similar way to characterise the frequency response of the beta of a rolling regression as a function of its window length in a similar manner?

Note that I have read this paper which explores using regression beta as a filter but I was unable to get an intuitive understanding of why the larger the window length is, the lower the cutoff frequency appears to be.

Here are a few descriptive charts of our data that might show why we are asking this question:

In the chart above, the thick grey line is the regression over all the data (i.e. a long window length) and has a slope of 1.79. The coloured lines (and points) show the regressions for single years of data (the betas for these lines are in the brackets in the legend of the chart). The beta for any individual year is significantly lower than the overall beta. In other words the grey line is much steeper than the other lines. This is true for any arbitrary 1 year window on the data (i.e. there is no year-long window of data that has a slope as steep as the overall data).

This next chart shows the rolling betas for different window lengths:

So the question we are trying to answer (and hoping to do so by characterising regression as a LPF) is why do the lines in this chart have higher biases (higher on average) as the window length increases (e.g. the yellow line is higher than the green line). How can we interpret this phenomenon?


texts: Estimate a microphone's transfer function

texts:
I am trying to compare the transfer functions of k different microphones and need some advice on my current approach. As of now, I have attempted the following:

Play a pink noise test tone $(X)$ from a speaker and record the tone on all $k$ target microphones. Let's call the outputs as $Y_{i}$.
Compute the inverse pink noise tone and convolve it with the recorded microphone tones $Y_{i}$ separately. The output of this convolution operation should be the estimated transfer function ($H_{i}$) of each microphone (it will not be 100% accurate because there is an effect of the room conditions, but let's ignore it for now).

Now to validate that the estimated transfer functions of all microphones are correct, I plan to compute the inverse transfer function $(H^{-1})$, multiply it with each microphone output $Y_i$ in the frequency domain and recover the input signal X. Finally, if all goes well, I expect the recovered X to be the same for all microphones.

$$ Y_i = H_i . X  $$

$$X = H_i^{-1} . Y_i $$

Are there any other approaches to do this? Any pointers for Python implementations of it?


texts: How Matched filter is transformed due to frequency modulation?

texts:
Intro:

Consider a signal source S which emits a PAM modulated signal, bits encoded as +g(t) and -g(t) pulses.

The matched filter for S is f - f is a function of g.

Then S is fed to an FM modulator and yields a modulated signal, S', let's call the corresponding matched filter of the modulated signal f'.

Question:

If f is known how can f' be computed ?


texts: Understanding a lowpass - comb filter implementation

texts:
I try to understand the implementation of the low-pass comb filter of the Freeverb reverberation algorithm:

https://ccrma.stanford.edu/~jos/pasp/Lowpass_Feedback_Comb_Filter.html

The original implementation is as follows:

inline float comb::process(float input)
{
    float output;
    output = buffer[bufidx];
    filterstore = (output*damp2) + (filterstore*damp1);
    buffer[bufidx] = input + (filterstore*feedback);
    if(++bufidx>=bufsize) bufidx = 0;
    return output;
}

I drew the closed loop of this algorithm:

In the link, it says:

  Inspection of comb.h in the Freeverb source shows that Freeverb's
  ``comb'' filter is more specifically a lowpass-feedback-comb filter
  (LBCF4.11--§2.6.5). It is constructed using a delay line whose output
  is lowpass-filtered and summed with the delay-line's input. The
  particular lowpass used in Freeverb is a unity-gain one-pole lowpass
  having the transfer function

  $\displaystyle H(z) = \frac{1-d}{1-d\,z^{-1}}. $

  When $ d=0$ , the LBCF reduces to the feedback comb filter (FBCF) of
  §2.6.2 in which the feedback was not filtered. The overall LBCF
  transfer function is then

  $ \displaystyle \hbox{LBCF}_{N}^{\,f,\,d} \;= \; \frac{z^{-N}}{1 -f\frac{1-d}{1-d\,z^{-1}}\,z^{-N}}. $

Apparently, this transfer function is implemented here. But can someone tell me, how to derive this implementation out of this transfer function?


texts:
I'm obtaining MFCC by using the function wave2mfcc from the Speech and Audio Processing library (SAP) from here http://mirlab.org/jang/matlab/toolbox/

The steps follows the procedure:

Take the Fourier transform of (a windowed excerpt of) a signal.
Map the powers of the spectrum obtained above onto the mel scale, using triangular overlapping windows.
Take the logs of the powers at each of the mel frequencies
Take the discrete cosine transform of the list of mel log powers, as if it were a signal

As I have previously read, the MFCC are the spectral envelope of the signal, but I don't know how to relate the coefficients with the envelope.

For example, if I have a signal 'example.wav'

[signal fs] = wavread('example.wav');
    figure(1); plot(log(abs(fft(signal)));
    mfcc = wave3mfcc(signal,fs);
    figure(2); plot(mfcc):

As I understand the MFCC should have the shape of the figure 1, but in my case it doesn't, maybe I'm missing some step.

Can anybody tell me please how to relate them and being able to graph?


texts: Is there a repository for source separation test images e.g. the star-galaxy images?

texts:
I'm having trouble finding standard images used for testing source-separation algorithms in image processing. For instance, a common example I see is the "star-galaxy" images. Most papers I've seen cite the following paper as their source of the image:

M. B. McCoy, V. Cevher, Q. T. Dinh, A. Asaei, and L. Baldassarre,Convexity in source separation:  Models, geometry, and algorithms, IEEE Signal Processing Magazine, 31 (2014), pp. 87–95.

However, upon looking at that article, it says their image is from NASA/JPL-Caltech and used with their permission. I haven't been able to find a download link anywhere. Any references in the right direction would be very helpful!


texts: Is a pulse of white noise still properly described as stationary?

texts:
I understand that a signal that is white noise is stationary (or more properly that the process generating it is stationary).

What if the white noise is delivered as a single pulse or a series of pulses with silence between them? Is that still stationary, or does the pulsatility suddenly convert this into a non-stationary system?

If it is no longer stationary, is there a term of art conveying the information that the signal within the pulse is stationary (e.g., a pure tone or white noise) rather than  a non-stationary natural or musical sound?